
import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt
import numpy as np
import mapclassify


def wf(x):
    return "/Users/andywang/desktop/Lab3/" + x
shp = gpd.read_file(wf("pc_MA_EE508.shp"))
d = pd.read_csv(wf("pc_MA_EE508.csv"))
shp_ma_outline = gpd.read_file(wf("MA_outline.shp"))

/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/3228652068.py:4: DtypeWarning: Columns (20) have mixed types. Specify dtype option on import or set low_memory=False.
  d = pd.read_csv(wf("pc_MA_EE508.csv"))


shp = shp.set_index('pid')
d = d.set_index('pid')


shp_c = shp.copy()
shp_c['geometry'] = shp_c['geometry'].centroid


shp_c['ms'] = d['ha'] / 10
shp_c.join(d['travel']).plot(column='travel', markersize='ms', alpha=0.5, scheme = 'quantiles')

<Axes: >


# Compute the change in percentage of parcel-level forest cover between 1990 and 2010
d['p_f_change'] = d['p_f_2010'] - d['p_f_1990']

# Compute the log of the per-hectare sale price of the last sale of the property
d['ls_price_per_ha_log'] = np.log(d['ls_price'] / d['ha'])

# Convert the 'ls_date' column to a datetime format and extract the year
d['ls_year'] = pd.to_datetime(d['ls_date']).dt.year

# Create a new variable 'p_e_ls' to represent the percentage of easement at the time of the last sale
# If 'e_year' is less than or equal to 'ls_year', use 'p_e'; otherwise, set it to 0
d['p_e_ls'] = np.where(d['e_year'] <= d['ls_year'], d['p_e'], 0)

# Normalize the building area to the parcel size
d['bld_area_per_ha'] = d['bld_area'] / d['ha']


# Define the list of predictor variables
XCOLS_LS = ['slope', 'travel', 'p_wet', 'coast_2500', 'river_frontage', 
            'lake_frontage', 'pop_dens_tract_2012-2016', 
            'hh_inc_avg_tract_2012-2016',
            'p_f_2000', 'p_prot_2000_500', 
            'p_d_2000_200', 'p_d_2000_1000', 'p_d_2000_5000',
            'p_dev_open', 'p_dev_low', 'p_dev_medium', 'p_dev_high', 
            'p_crops', 'p_pasture','p_grassland', 'bld_area_per_ha', 
            'ls_year']

# Create the X and y DataFrames/Series for last sale price prediction
X_ls = d[d['ls_year'] >= 2000][XCOLS_LS]  # Select data from 2000 and onward
y_ls = d[d['ls_year'] >= 2000]['ls_price_per_ha_log']  # Sale price per hectare with log transformation


# Define the list of predictor variables
XCOLS_FC = ['slope', 'travel', 'p_wet', 'coast_2500', 'river_frontage', 
            'lake_frontage',  'pop_dens_tract_2012-2016', 
            'hh_inc_avg_tract_2012-2016', 
            'p_f_1990', 'p_prot_1990_500', 
            'p_d_1990_200', 'p_d_1990_1000', 'p_d_1990_5000']

# Create the X and y DataFrames/Series for forest change probability estimation
X_fc = d[XCOLS_FC]  # All data without any time filter
y_fc = d['p_f_change']  # Change in percentage of forest cover from 1990 to 2010


# Define the selection criteria for X_ls and y_ls
i_ls = d[XCOLS_LS].isnull().sum(1).eq(0)
i_ls &= d['ls_year'].ge(2000)
i_ls &= d['ls_year'].le(2019)  # Sold between 2000 and 2019
i_ls &= d['p_e_ls'].eq(0)  # Not protected by an easement at the time of sale

# Select rows based on the criteria
X_ls = d.loc[i_ls, XCOLS_LS].copy()
y_ls = d.loc[i_ls, 'ls_price_per_ha_log'].copy()

# Convert 'ls_year' values to integers
X_ls['ls_year'] = X_ls['ls_year'].astype(int)


# Define the selection criteria for X_fc and y_fc
i_fc = d[XCOLS_FC].isnull().sum(1).eq(0)
i_fc &= d['p_prot'].le(5)  # Not protected (p_prot <= 5%)

# Select rows based on the criteria
X_fc = d.loc[i_fc, XCOLS_FC].copy()
y_fc = d.loc[i_fc, 'p_f_change'].copy()


import statsmodels.api as sm
model_ls_ols = sm.OLS(y_ls, sm.add_constant(X_ls)).fit()
model_ls_ols.summary()


fig, ax = plt.subplots(figsize=(17, 10))
dp = shp_c.join(model_ls_ols.fittedvalues.rename('pred'), how='inner')
shp_ma_outline.plot(ax=ax, color='none', edgecolor='black')
dp.plot(
    'pred',
    ax=ax,
    markersize='ms',
    alpha=0.5,
    cmap='Spectral_r',
    scheme='quantiles',
    k=10,
    legend=True
)
ax.axis('off')

(19014.26383461053, 345702.9467683639, 768457.3790219357, 968851.8563356873)


fig, ax = plt.subplots(figsize=(17, 10))
rp = shp_c.join(model_ls_ols.resid.rename('resid'), how='inner')
shp_ma_outline.plot(ax=ax, color='none', edgecolor='black')
rp.plot(
    'resid',
    ax=ax,
    markersize='ms',
    alpha=0.5,
    cmap='Spectral_r',
    scheme='quantiles',
    k=10,
    legend=True
)
ax.axis('off')

(19014.26383461053, 345702.9467683639, 768457.3790219357, 968851.8563356873)


df = X_ls[['ls_year']].join(model_ls_ols.resid.rename('resid'))
df.groupby('ls_year')['resid'].mean().plot()

<Axes: xlabel='ls_year'>


from sklearn.metrics import mean_squared_error
mean_squared_error(y_ls, model_ls_ols.fittedvalues)

1.980329829753318


from sklearn.preprocessing import OneHotEncoder

enc = OneHotEncoder(categories='auto').fit(X_ls[['ls_year']])
dummies_array = enc.transform(X_ls[['ls_year']]).toarray()
dummies_df = pd.DataFrame(
    dummies_array,
    index=X_ls.index,
    columns=[v[8:] for v in enc.get_feature_names_out()],
)


X_ls2 = X_ls.drop(columns='ls_year').join(dummies_df).drop(columns='2000')


model_ls_ols2 = sm.OLS(y_ls, sm.add_constant(X_ls2)).fit()
model_ls_ols2.summary()


mean_squared_error(y_ls, model_ls_ols2.fittedvalues)

1.9570070916504736


model_ls_ols2.params.loc['2001':'2017'].plot()

<Axes: >


from itertools import combinations
X_ls3 = X_ls.drop(columns='ls_year').copy()
cols = X_ls3.columns
for v in cols:
    X_ls3[v + '^2'] = X_ls3[v] ** 2
for v1, v2 in combinations(cols,2):
    X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
X_ls3 = X_ls3.join(dummies_df)

/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]
/var/folders/tk/v2hjsgjj3ysfqpw72wqmt1m40000gn/T/ipykernel_14366/1670456788.py:7: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  X_ls3[v1 + '*' + v2] = X_ls3[v1] * X_ls3[v2]


model_ls_ols3 = sm.OLS(y_ls, sm.add_constant(X_ls3)).fit()
model_ls_ols3.summary()


mean_squared_error(y_ls, model_ls_ols3.fittedvalues)

1.6841173295590273


fig, ax = plt.subplots(figsize=(17, 10))
rp3 = shp_c.join(model_ls_ols3.resid.rename('resid'), how='inner')
shp_ma_outline.plot(ax=ax, color='none', edgecolor='black')
rp3.plot(
    'resid',
    ax=ax,
    markersize='ms',
    alpha=0.5,
    cmap='Spectral_r',
    scheme='quantiles',
    k=10,
    legend=True
)
ax.set_title("Residuals of Land Sale Prices Estimation in Massachusetts", fontsize=18)

ax.axis('off')
#plt.savefig("/Users/andywang/desktop/1_ls_ols3_resid.png", dpi=200)

(19014.26383461053, 345702.9467683639, 768457.3790219357, 968851.8563356873)


fig, ax = plt.subplots(figsize=(17, 10))
dp3 = shp_c.join(model_ls_ols3.fittedvalues.rename('pred'), how='inner')
shp_ma_outline.plot(ax=ax, color='none', edgecolor='black')
dp3.plot(
    'pred',
    ax=ax,
    markersize='ms',
    alpha=0.5,
    cmap='Spectral_r',
    scheme='quantiles',
    k=10,
    legend=True
)
ax.set_title("Predicted Values of Land Sale Prices Estimation in Massachusetts", fontsize=18)
ax.axis('off')
#plt.savefig("/Users/andywang/desktop/1_ls_ols3_pred.png", dpi=200)

(19014.26383461053, 345702.9467683639, 768457.3790219357, 968851.8563356873)


from sklearn.linear_model import LinearRegression
from sklearn.ensemble import (
    ExtraTreesRegressor, 
    GradientBoostingRegressor, 
    RandomForestRegressor
)
from sklearn.model_selection import train_test_split, KFold


model = LinearRegression()
model.fit(X_ls, y_ls)

LinearRegression()

LinearRegression()


model.predict(X_ls)
mean_squared_error(y_ls, model.predict(X_ls))

1.980329829753318


X_train, X_test, y_train, y_test = train_test_split(X_ls, y_ls, test_size=0.5)
model = LinearRegression()
model.fit(X_train, y_train)
mean_squared_error(y_test, model.predict(X_test))

2.011834891126498


X = X_ls
y = y_ls
mse_train = []
mse_test = []
for i_train, i_test in KFold(n_splits=5, shuffle=True).split(X):
    X_train, X_test = X.iloc[i_train], X.iloc[i_test]
    y_train, y_test = y.iloc[i_train], y.iloc[i_test]
    model = LinearRegression()
    model.fit(X_train, y_train)
    mse_train += [mean_squared_error(y_train, model.predict(X_train))]
    mse_test += [mean_squared_error(y_test, model.predict(X_test))]
print('Avg MSE (within-sample):', round(np.mean(mse_train), 4))
print('Avg MSE (out-of-sample):', round(np.mean(mse_test), 4))

Avg MSE (within-sample): 1.9799
Avg MSE (out-of-sample): 1.984


X = X_ls2
y = y_ls
mse_train = []
mse_test = []
for i_train, i_test in KFold(n_splits=5, shuffle=True).split(X):
    X_train, X_test = X.iloc[i_train], X.iloc[i_test]
    y_train, y_test = y.iloc[i_train], y.iloc[i_test]
    model = LinearRegression()
    model.fit(X_train, y_train)
    mse_train += [mean_squared_error(y_train, model.predict(X_train))]
    mse_test += [mean_squared_error(y_test, model.predict(X_test))]
print('Avg MSE (within-sample):', round(np.mean(mse_train), 4))
print('Avg MSE (out-of-sample):', round(np.mean(mse_test), 4))

Avg MSE (within-sample): 1.9565
Avg MSE (out-of-sample): 1.9621


X = X_ls3
y = y_ls
mse_train = []
mse_test = []
for i_train, i_test in KFold(n_splits=5, shuffle=True).split(X):
    X_train, X_test = X.iloc[i_train], X.iloc[i_test]
    y_train, y_test = y.iloc[i_train], y.iloc[i_test]
    model = LinearRegression()
    model.fit(X_train, y_train)
    mse_train += [mean_squared_error(y_train, model.predict(X_train))]
    mse_test += [mean_squared_error(y_test, model.predict(X_test))]
print('Avg MSE (within-sample):', round(np.mean(mse_train), 4))
print('Avg MSE (out-of-sample):', round(np.mean(mse_test), 4))

Avg MSE (within-sample): 1.6781
Avg MSE (out-of-sample): 1.7686


X = X_ls
y = y_ls
mse_train = []
mse_test = []
for i_train, i_test in KFold(n_splits=5, shuffle=True).split(X):
    X_train, X_test = X.iloc[i_train], X.iloc[i_test]
    y_train, y_test = y.iloc[i_train], y.iloc[i_test]
    model = GradientBoostingRegressor()
    model.fit(X_train, y_train)
    mse_train += [mean_squared_error(y_train, model.predict(X_train))]
    mse_test += [mean_squared_error(y_test, model.predict(X_test))]
print('Avg MSE (within-sample):', round(np.mean(mse_train), 4))
print('Avg MSE (out-of-sample):', round(np.mean(mse_test), 4))

Avg MSE (within-sample): 1.4143
Avg MSE (out-of-sample): 1.4732


X = X_ls
y = y_ls
mse_train = []
mse_test = []
for i_train, i_test in KFold(n_splits=5, shuffle=True).split(X):
    X_train, X_test = X.iloc[i_train], X.iloc[i_test]
    y_train, y_test = y.iloc[i_train], y.iloc[i_test]
    model = RandomForestRegressor(n_estimators=20)
    model.fit(X_train, y_train)
    mse_train += [mean_squared_error(y_train, model.predict(X_train))]
    mse_test += [mean_squared_error(y_test, model.predict(X_test))]
print('Avg MSE (within-sample):', round(np.mean(mse_train), 4))
print('Avg MSE (out-of-sample):', round(np.mean(mse_test), 4))

Avg MSE (within-sample): 0.2256
Avg MSE (out-of-sample): 1.4086


X = X_ls
y = y_ls
mse_train = []
mse_test = []
for i_train, i_test in KFold(n_splits=5, shuffle=True).split(X):
    X_train, X_test = X.iloc[i_train], X.iloc[i_test]
    y_train, y_test = y.iloc[i_train], y.iloc[i_test]
    model = ExtraTreesRegressor(n_estimators=50)
    model.fit(X_train, y_train)
    mse_train += [mean_squared_error(y_train, model.predict(X_train))]
    mse_test += [mean_squared_error(y_test, model.predict(X_test))]
print('Avg MSE (within-sample):', round(np.mean(mse_train), 4))
print('Avg MSE (out-of-sample):', round(np.mean(mse_test), 4))

Avg MSE (within-sample): 0.0
Avg MSE (out-of-sample): 1.3451


feat_imp = pd.Series(model.feature_importances_, index=X.columns)
feat_imp.sort_values().plot(kind='barh')

<Axes: >


fig, ax = plt.subplots(figsize=(17, 10))
y_pred = pd.Series(model.predict(X), index=X.index)
resid = (y - y_pred).rename('resid')
shp_ma_outline.plot(ax=ax, color='none', edgecolor='black')
shp_p = shp_c.join(resid, how='inner')
shp_p.plot(
    'resid',
    ax=ax,
    markersize='ms',
    scheme='quantiles',
    legend=True,
    cmap='Spectral',
    figsize=(17, 10)
)

<Axes: >


model = ExtraTreesRegressor(n_estimators=50, bootstrap=True, oob_score=True)
model.fit(X, y)
y_pred = pd.Series(model.oob_prediction_, index=X.index)
resid = (y - y_pred).rename('resid')
shp_p = shp_c.join(resid, how='inner')
fig, ax = plt.subplots(figsize=(17, 10))
shp_ma_outline.plot(ax=ax, color='none', edgecolor='black')
shp_p.plot(
    'resid',
    ax = ax,
    markersize='ms',
    scheme='quantiles',
    legend=True,
    cmap='Spectral',
    figsize=(17, 10)
)

<Axes: >


model = RandomForestRegressor(n_estimators=50, bootstrap=True, oob_score=True)
model.fit(X, y)
y_pred = pd.Series(model.oob_prediction_, index=X.index)
resid = (y - y_pred).rename('resid')
shp_p = shp_c.join(resid, how='inner')
fig, ax = plt.subplots(figsize=(17, 10))
shp_ma_outline.plot(ax=ax, color='none', edgecolor='black')
shp_p.plot(
    'resid',
    ax = ax,
    markersize='ms',
    scheme='quantiles',
    legend=True,
    cmap='Spectral',
    figsize=(17, 10)
)

<Axes: >


shp_c['x'] = shp_c['geometry'].x
shp_c['y'] = shp_c['geometry'].y
coords = shp_c[['x', 'y']]


X = X_ls.join(coords)
y = y_ls
mse_train = []
mse_test = []
for i_train, i_test in KFold(n_splits=5, shuffle=True).split(X):
    X_train, X_test = X.iloc[i_train], X.iloc[i_test]
    y_train, y_test = y.iloc[i_train], y.iloc[i_test]
    model = GradientBoostingRegressor()
    model.fit(X_train, y_train)
    mse_train += [mean_squared_error(y_train, model.predict(X_train))]
    mse_test += [mean_squared_error(y_test, model.predict(X_test))]
print('Avg MSE (within-sample):', round(np.mean(mse_train), 4))
print('Avg MSE (out-of-sample):', round(np.mean(mse_test), 4))

Avg MSE (within-sample): 1.37
Avg MSE (out-of-sample): 1.4342


X = X_ls.join(coords)
y = y_ls
mse_train = []
mse_test = []
for i_train, i_test in KFold(n_splits=5, shuffle=True).split(X):
    X_train, X_test = X.iloc[i_train], X.iloc[i_test]
    y_train, y_test = y.iloc[i_train], y.iloc[i_test]
    model = RandomForestRegressor(n_estimators=20)
    model.fit(X_train, y_train)
    mse_train += [mean_squared_error(y_train, model.predict(X_train))]
    mse_test += [mean_squared_error(y_test, model.predict(X_test))]
print('Avg MSE (within-sample):', round(np.mean(mse_train), 4))
print('Avg MSE (out-of-sample):', round(np.mean(mse_test), 4))

Avg MSE (within-sample): 0.2167
Avg MSE (out-of-sample): 1.3695


X = X_ls.join(coords)
y = y_ls
mse_train = []
mse_test = []
for i_train, i_test in KFold(n_splits=5, shuffle=True).split(X):
    X_train, X_test = X.iloc[i_train], X.iloc[i_test]
    y_train, y_test = y.iloc[i_train], y.iloc[i_test]
    model = ExtraTreesRegressor(n_estimators=50)
    model.fit(X_train, y_train)
    mse_train += [mean_squared_error(y_train, model.predict(X_train))]
    mse_test += [mean_squared_error(y_test, model.predict(X_test))]
print('Avg MSE (within-sample):', round(np.mean(mse_train), 4))
print('Avg MSE (out-of-sample):', round(np.mean(mse_test), 4))

Avg MSE (within-sample): 0.0
Avg MSE (out-of-sample): 1.283


X = X_ls.join(coords)
y = y_ls
mse_train = []
mse_test = []
for i_train, i_test in KFold(n_splits=5, shuffle=True).split(X):
    X_train, X_test = X.iloc[i_train], X.iloc[i_test]
    y_train, y_test = y.iloc[i_train], y.iloc[i_test]
    model = GradientBoostingRegressor(n_estimators=1000)
    model.fit(X_train, y_train)
    mse_train += [mean_squared_error(y_train, model.predict(X_train))]
    mse_test += [mean_squared_error(y_test, model.predict(X_test))]
print('Avg MSE (within-sample):', round(np.mean(mse_train), 4))
print('Avg MSE (out-of-sample):', round(np.mean(mse_test), 4))

Avg MSE (within-sample): 1.0249
Avg MSE (out-of-sample): 1.3143


X = X_ls.join(coords)
y = y_ls
mse_train = []
mse_test = []
for i_train, i_test in KFold(n_splits=5, shuffle=True).split(X):
    X_train, X_test = X.iloc[i_train], X.iloc[i_test]
    y_train, y_test = y.iloc[i_train], y.iloc[i_test]
    model = RandomForestRegressor(n_estimators=100, oob_score=True)
    model.fit(X_train, y_train)
    mse_train += [mean_squared_error(y_train, model.predict(X_train))]
    mse_test += [mean_squared_error(y_test, model.predict(X_test))]
print('Avg MSE (within-sample):', round(np.mean(mse_train), 4))
print('Avg MSE (out-of-sample):', round(np.mean(mse_test), 4))

Avg MSE (within-sample): 0.1833
Avg MSE (out-of-sample): 1.3001


X = X_ls.join(coords)
y = y_ls
mse_train = []
mse_test = []
for i_train, i_test in KFold(n_splits=5, shuffle=True).split(X):
    X_train, X_test = X.iloc[i_train], X.iloc[i_test]
    y_train, y_test = y.iloc[i_train], y.iloc[i_test]
    model = ExtraTreesRegressor(
        n_estimators=100, bootstrap=True, oob_score=True
    )
    model.fit(X_train, y_train)
    mse_train += [mean_squared_error(y_train, model.predict(X_train))]
    mse_test += [mean_squared_error(y_test, model.predict(X_test))]
print('Avg MSE (within-sample):', round(np.mean(mse_train), 4))
print('Avg MSE (out-of-sample):', round(np.mean(mse_test), 4))

Avg MSE (within-sample): 0.1819
Avg MSE (out-of-sample): 1.2904


X = X_ls.join(coords)
y = y_ls
mse_train = []
mse_test = []
for i_train, i_test in KFold(n_splits=5, shuffle=True).split(X):
    X_train, X_test = X.iloc[i_train], X.iloc[i_test]
    y_train, y_test = y.iloc[i_train], y.iloc[i_test]
    model = ExtraTreesRegressor(
        n_estimators=250, bootstrap=True, oob_score=True, n_jobs=-1
    )
    model.fit(X_train, y_train)
    mse_train += [mean_squared_error(y_train, model.predict(X_train))]
    mse_test += [mean_squared_error(y_test, model.predict(X_test))]
print('Avg MSE (within-sample):', round(np.mean(mse_train), 4))
print('Avg MSE (out-of-sample):', round(np.mean(mse_test), 4))

Avg MSE (within-sample): 0.1767
Avg MSE (out-of-sample): 1.2843


xcols_ls_pred = [v for v in XCOLS_LS if v != 'ls_year']
i_ls_pred = d[xcols_ls_pred].isnull().sum(1).eq(0) & d['p_prot'].le(5)
X_ls_pred = d.loc[i_ls_pred, xcols_ls_pred]
X_ls_pred['ls_year'] = 2010


model_et_coords_250 = ExtraTreesRegressor(
    n_estimators=250, n_jobs=-1, bootstrap=True, oob_score=True
) 
model_et_coords_250.fit(X_ls.join(coords), y_ls)
y_pred_et_coords_250 = pd.Series(
    model_et_coords_250.predict(X_ls_pred.join(coords)), index=X_ls_pred.index
)


fig, ax = plt.subplots(figsize=(17, 10))
dp_et = shp_c.join(y_pred_et_coords_250.rename('pred1'), how='inner')
shp_ma_outline.plot(ax=ax, color='none', edgecolor='black')
dp_et.plot(
    'pred1',
    ax=ax,
    markersize='ms',
    alpha=0.5,
    cmap='Spectral_r',
    scheme='quantiles',
    k=10,
    legend=True
)
ax.set_title("Predicted Values of Land Sale Prices Estimation in Massachusetts", fontsize=18)
ax.axis('off')
#plt.savefig("/Users/andywang/desktop/1_ls_et_coords_250_pred.png", dpi=200)

(19014.26383461053, 345702.9467683639, 768457.3790219357, 968851.8563356873)


resid = (y - y_pred_et_coords_250).rename('resid')
shp_p = shp_c.join(resid, how='inner')
fig, ax = plt.subplots(figsize=(17, 10))
shp_ma_outline.plot(ax=ax, color='none', edgecolor='black')
shp_p.plot(
    'resid',
    ax = ax,
    markersize='ms',
    scheme='quantiles',
    legend=True,
    cmap='Spectral',
    figsize=(17, 10)
)
ax.set_title("Residuals of Land Sale Prices Estimation in Massachusetts", fontsize=18)
ax.axis('off')
#plt.savefig("/Users/andywang/desktop/1_ls_et_coords_250_resid.png", dpi=200)

(19014.26383461053, 345702.9467683639, 768457.3790219357, 968851.8563356873)


model_rf_coords_100 = RandomForestRegressor(
    n_estimators=100, n_jobs=-1, oob_score=True
) 
model_rf_coords_100.fit(X_ls.join(coords), y_ls)
y_pred_rf_coords_100 = pd.Series(
    model_rf_coords_100.predict(X_ls_pred.join(coords)), index=X_ls_pred.index
)


fig, ax = plt.subplots(figsize=(17, 10))
dp_rf = shp_c.join(y_pred_rf_coords_100.rename('pred2'), how='inner')
shp_ma_outline.plot(ax=ax, color='none', edgecolor='black')
dp_rf.plot(
    'pred2',
    ax=ax,
    markersize='ms',
    alpha=0.5,
    cmap='Spectral_r',
    scheme='quantiles',
    k=10,
    legend=True
)
ax.set_title("Predicted Values of Land Sale Prices Estimation in Massachusetts", fontsize=18)
ax.axis('off')
#plt.savefig("/Users/andywang/desktop/1_ls_rf_coords_100_pred.png", dpi=200)

(19014.26383461053, 345702.9467683639, 768457.3790219357, 968851.8563356873)


resid = (y - y_pred_rf_coords_100).rename('resid')
shp_p = shp_c.join(resid, how='inner')
fig, ax = plt.subplots(figsize=(17, 10))
shp_ma_outline.plot(ax=ax, color='none', edgecolor='black')
shp_p.plot(
    'resid',
    ax = ax,
    markersize='ms',
    scheme='quantiles',
    legend=True,
    cmap='Spectral',
    figsize=(17, 10)
)
ax.set_title("Residuals of Land Sale Prices Estimation in Massachusetts", fontsize=18)
ax.axis('off')
#plt.savefig("/Users/andywang/desktop/1_ls_rf_coords_100_resid.png", dpi=200)

(19014.26383461053, 345702.9467683639, 768457.3790219357, 968851.8563356873)


xcols_fc_pred = [v for v in XCOLS_FC if v != 'fc_year']
i_fc_pred = d[xcols_fc_pred].isnull().sum(1).eq(0) & d['p_prot'].le(5)
X_fc_pred = d.loc[i_fc_pred, xcols_fc_pred]
model_fc_coords_250 = ExtraTreesRegressor(
    n_estimators=250, n_jobs=-1, bootstrap=True, oob_score=True
)
model_fc_coords_250.fit(X_fc.join(coords), y_fc)
y_pred_fc_coords_250 = pd.Series(
    model_fc_coords_250.predict(X_fc_pred.join(coords)), index=X_fc_pred.index
)


fig, ax = plt.subplots(figsize=(17, 10))
dp_fc = shp_c.join(y_pred_fc_coords_250.rename('pred3'), how='inner')

shp_ma_outline.plot(ax=ax, color='none', edgecolor='black')
dp_fc.plot(
    'pred3',
    ax=ax,
    markersize='ms',
    alpha=0.5,
    cmap='Spectral', 
    scheme='quantiles',
    k=10,
    legend=True
)
ax.set_title("Predicted Forest Cover Change in Massachusetts", fontsize=18)
ax.axis('off')
#plt.savefig("/Users/andywang/desktop/1_fc_et_pred.png", dpi=200)

(19014.26383461053, 345702.9467683639, 768457.3790219357, 968851.8563356873)


yp = pd.read_csv(wf("y_predictions.csv"))
yp = yp.set_index('pid')


dm = d.join(yp, how='inner')
dm


dm['p_f_loss_pred_100yr'] = dm['fc_pred'].mul(-10).clip(upper=dm['p_f_2010'], lower=0)
dm['ha_f_loss_pred_100yr'] = dm['p_f_loss_pred_100yr'] * dm['ha'] / 100
dm['tco2_loss_pred_100yr'] = dm['ha_f_loss_pred_100yr'] * 400
dm['acq_cost_pred_usd'] = np.exp(dm['ls_pred']) * dm['ha']
dm['ce_cost_pred_usd'] = dm['acq_cost_pred_usd'] * 0.4
dm['ce_cost_$/tco2'] = dm['ce_cost_pred_usd'] / dm['tco2_loss_pred_100yr']


dm_filtered = dm[~dm['ce_cost_$/tco2'].isin([np.inf, -np.inf])]


fig, ax = plt.subplots(figsize=(17, 10))
shp_ma_outline.plot(ax=ax, color='none', edgecolor='black')
dm0 = shp_c.join(dm_filtered['ce_cost_$/tco2'].rename('avoidedCO2'), how='inner')
dm0.plot(
    'avoidedCO2',
    ax=ax,
    markersize='ms',
    alpha=0.5,
    cmap='Spectral',
    scheme='quantiles',
    k=10,
    legend=True
)
ax.set_title("Conservation Easement Cost per Ton of Expected Avoided CO2 Emissions", fontsize=18)
ax.axis('off')
#plt.savefig("/Users/andywang/desktop/2_ce_cost_usd_per_tco2.png", dpi=200)

(19014.26383461053, 345702.9467683639, 768457.3790219357, 968851.8563356873)


dm_sorted = dm_filtered.sort_values(by='ce_cost_$/tco2')
dm_sorted['cum_tco2'] = dm_sorted['tco2_loss_pred_100yr'].cumsum()
dm_supply_curve = dm_sorted[dm_sorted['ce_cost_$/tco2'] <= 1000]
dm_sorted


fig, ax = plt.subplots(figsize=(10, 6))
dm_supply_curve.plot('cum_tco2', 'ce_cost_$/tco2', ax=ax)
ax.set_title("Estimated Supply Curve for Avoided CO2 Emissions", fontsize=18)
ax.set_xlabel("Total Avoided CO2 Emissions (tons)")
ax.set_ylabel("Cost per Ton of Avoided CO2 ($/tCO2)")
#plt.savefig("/Users/andywang/desktop/2_tco2_supply_curve.png", dpi=200)

Text(0, 0.5, 'Cost per Ton of Avoided CO2 ($/tCO2)')


rggi_price = 13
scc_price = 51
uruguay_price = 135
def outcome_variables(dm, carbon_price):
    selected_parcels = dm[dm['ce_cost_$/tco2'] <= carbon_price]
    n = len(selected_parcels)
    ha_pc = selected_parcels['ha'].sum()
    ha_loss_avo = selected_parcels['ha_f_loss_pred_100yr'].sum()
    tco2_avo = selected_parcels['tco2_loss_pred_100yr'].sum()
    budget = selected_parcels['ce_cost_pred_usd'].sum()
    avg_usd_tco2 = budget / tco2_avo
    return n, ha_pc, ha_loss_avo, tco2_avo, budget, avg_usd_tco2
rggi_outcomes = outcome_variables(dm, rggi_price)
scc_outcomes = outcome_variables(dm, scc_price)
uruguay_outcomes = outcome_variables(dm, uruguay_price)


scenario_df = pd.DataFrame({
    'Scenario': [rggi_price, scc_price, uruguay_price],
    'n': [rggi_outcomes[0], scc_outcomes[0], uruguay_outcomes[0]],
    'ha_pc': [rggi_outcomes[1], scc_outcomes[1], uruguay_outcomes[1]],
    'ha_loss_avo': [rggi_outcomes[2], scc_outcomes[2], uruguay_outcomes[2]],
    'tco2_avo': [rggi_outcomes[3], scc_outcomes[3], uruguay_outcomes[3]],
    'budget': [rggi_outcomes[4], scc_outcomes[4], uruguay_outcomes[4]],
    'avg_usd_tco2': [rggi_outcomes[5], scc_outcomes[5], uruguay_outcomes[5]]
})
#scenario_df.to_csv("/Users/andywang/desktop/2_carbon_prices_summary.csv", index=False)


BUDGET=70000000
dm_sorted = dm.sort_values(by='ce_cost_$/tco2')
dm_sorted['ce_cost_pred_usd_sum'] = dm_sorted['ce_cost_pred_usd'].cumsum()
i_scenario_1 = dm_sorted['ce_cost_pred_usd_sum'].lt(BUDGET)
selected_parcels_scenario_1 = dm_sorted[i_scenario_1]


n_scenario_1 = len(selected_parcels_scenario_1)
ha_pc_scenario_1 = selected_parcels_scenario_1['ha'].sum()
ha_loss_avo_scenario_1 = selected_parcels_scenario_1['ha_f_loss_pred_100yr'].sum()
tco2_avo_scenario_1 = selected_parcels_scenario_1['tco2_loss_pred_100yr'].sum()
budget_scenario_1 = selected_parcels_scenario_1['ce_cost_pred_usd'].sum()
avg_usd_tco2_scenario_1 = budget_scenario_1 / tco2_avo_scenario_1


scenario_1_values = {
    'Scenario': 'Scenario 1',
    'n': n_scenario_1,
    'ha_pc': ha_pc_scenario_1,
    'ha_loss_avo': ha_loss_avo_scenario_1,
    'tco2_avo': tco2_avo_scenario_1,
    'budget': budget_scenario_1,
    'avg_usd_tco2': avg_usd_tco2_scenario_1
}
#pd.DataFrame([scenario_1_values]).to_csv("/Users/andywang/desktop/2_carbon_prices_summary.csv", index=False, mode='a', header=False)


fig1, ax1 = plt.subplots(figsize=(17, 10))
shp_ma_outline.plot(ax=ax1, color='none', edgecolor='black')
selected_parcels_scenario_1_gdf = shp_c.join(selected_parcels_scenario_1, on='pid')
selected_parcels_scenario_1_gdf.plot(ax=ax1, markersize='ha', color='blue', alpha=0.7)
ax1.axis('off')
ax1.set_title("Scenario 1", fontsize=18)
#plt.savefig("/Users/andywang/desktop/2_selected_parcels_scenario_1.png", dpi=200)

Text(0.5, 1.0, 'Scenario 1')


dm2010 = dm[dm['p_f_2010']>0]
dm_sorted2 = dm2010.sort_values(by='fc_pred')
dm_sorted2['ce_cost_pred_usd_sum'] = dm_sorted2['ce_cost_pred_usd'].cumsum()
i_scenario_2 = dm_sorted2['ce_cost_pred_usd_sum'].lt(BUDGET)
selected_parcels_scenario_2 = dm_sorted2[i_scenario_2]


n_scenario_2 = len(selected_parcels_scenario_2)
ha_pc_scenario_2 = selected_parcels_scenario_2['ha'].sum()
ha_loss_avo_scenario_2 = selected_parcels_scenario_2['ha_f_loss_pred_100yr'].sum()
tco2_avo_scenario_2 = selected_parcels_scenario_2['tco2_loss_pred_100yr'].sum()
budget_scenario_2 = selected_parcels_scenario_2['ce_cost_pred_usd'].sum()
avg_usd_tco2_scenario_2 = budget_scenario_2 / tco2_avo_scenario_2
scenario_2_values = {
    'Scenario': 'Scenario 2',
    'n': n_scenario_2,
    'ha_pc': ha_pc_scenario_2,
    'ha_loss_avo': ha_loss_avo_scenario_2,
    'tco2_avo': tco2_avo_scenario_2,
    'budget': budget_scenario_2,
    'avg_usd_tco2': avg_usd_tco2_scenario_2
}
#pd.DataFrame([scenario_2_values]).to_csv("/Users/andywang/desktop/2_carbon_prices_summary.csv", index=False, mode='a', header=False)


fig2, ax2 = plt.subplots(figsize=(17, 10))
shp_ma_outline.plot(ax=ax2, color='none', edgecolor='black')
selected_parcels_scenario_2_gdf = shp_c.join(selected_parcels_scenario_2, on='pid')
selected_parcels_scenario_2_gdf.plot(ax=ax2, markersize='ha', color='blue', alpha=0.7)
ax2.axis('off')
ax2.set_title("Scenario 2", fontsize=18)
#plt.savefig("/Users/andywang/desktop/2_selected_parcels_scenario_2.png", dpi=200)

Text(0.5, 1.0, 'Scenario 2')


dm['offer3'] = dm['p_f_2010']/100 * dm['ha'] * 2680
i_scenario_3 = dm['ce_cost_pred_usd'] < dm['offer3']
dm[i_scenario_3]['offer3'].sum()
selected_parcels_scenario_3 = dm[i_scenario_3]
dm['offer3']

pid
F_803174_3027860      1402.9264
M_238118_942044       7933.9792
M_240445_943056       3423.2980
M_245474_951109          0.0000
F_778757_3113679     18758.1240
                       ...     
M_74987_889210      232513.4232
F_182440_2984399     82739.6936
M_74001_894053      218594.6020
M_67459_886161      238781.4876
F_172467_2984485     70913.7380
Name: offer3, Length: 177277, dtype: float64


n_scenario_3 = len(selected_parcels_scenario_3)
ha_pc_scenario_3 = selected_parcels_scenario_3['ha'].sum()
ha_loss_avo_scenario_3 = selected_parcels_scenario_3['ha_f_loss_pred_100yr'].sum()
tco2_avo_scenario_3 = selected_parcels_scenario_3['tco2_loss_pred_100yr'].sum()
budget_scenario_3 = dm[i_scenario_3]['offer3'].sum()
avg_usd_tco2_scenario_3 = budget_scenario_3 / tco2_avo_scenario_3
scenario_3_values = {
    'Scenario': 'Scenario 3',
    'n': n_scenario_3,
    'ha_pc': ha_pc_scenario_3,
    'ha_loss_avo': ha_loss_avo_scenario_3,
    'tco2_avo': tco2_avo_scenario_3,
    'budget': budget_scenario_3,
    'avg_usd_tco2': avg_usd_tco2_scenario_3
}
#pd.DataFrame([scenario_3_values]).to_csv("/Users/andywang/desktop/2_carbon_prices_summary.csv", index=False, mode='a', header=False)


fig3, ax3 = plt.subplots(figsize=(17, 10))
shp_ma_outline.plot(ax=ax3, color='none', edgecolor='black')
selected_parcels_scenario_3_gdf = shp_c.join(selected_parcels_scenario_3, on='pid')
selected_parcels_scenario_3_gdf.plot(ax=ax3, markersize='ha', color='blue', alpha=0.7)
ax3.axis('off')
ax3.set_title("Scenario 3", fontsize=18)
#plt.savefig("/Users/andywang/desktop/2_selected_parcels_scenario_3.png", dpi=200)

Text(0.5, 1.0, 'Scenario 3')


dm['offer4']= dm['tco2_loss_pred_100yr'] *41
i_scenario_4 = dm['ce_cost_pred_usd'] < dm['offer4']
dm[i_scenario_4]['offer4'].sum()
selected_parcels_scenario_4 = dm[i_scenario_4]


n_scenario_4 = len(selected_parcels_scenario_4)
ha_pc_scenario_4 = selected_parcels_scenario_4['ha'].sum()
ha_loss_avo_scenario_4 = selected_parcels_scenario_4['ha_f_loss_pred_100yr'].sum()
tco2_avo_scenario_4 = selected_parcels_scenario_4['tco2_loss_pred_100yr'].sum()
budget_scenario_4 = dm[i_scenario_4]['offer4'].sum()
avg_usd_tco2_scenario_4 = budget_scenario_4 / tco2_avo_scenario_4
scenario_4_values = {
    'Scenario': 'Scenario 4',
    'n': n_scenario_4,
    'ha_pc': ha_pc_scenario_4,
    'ha_loss_avo': ha_loss_avo_scenario_4,
    'tco2_avo': tco2_avo_scenario_4,
    'budget': budget_scenario_4,
    'avg_usd_tco2': avg_usd_tco2_scenario_4
}
#pd.DataFrame([scenario_4_values]).to_csv("/Users/andywang/desktop/2_carbon_prices_summary.csv", index=False, mode='a', header=False)


fig4, ax4 = plt.subplots(figsize=(17, 10))
shp_ma_outline.plot(ax=ax4, color='none', edgecolor='black')
selected_parcels_scenario_4_gdf = shp_c.join(selected_parcels_scenario_4, on='pid')
selected_parcels_scenario_4_gdf.plot(ax=ax4, markersize='ha', color='blue', alpha=0.7)
ax4.axis('off')
ax4.set_title("Scenario 4", fontsize=18)
#plt.savefig("/Users/andywang/desktop/2_selected_parcels_scenario_4.png", dpi=200)

Text(0.5, 1.0, 'Scenario 4')


transaction_cost_per_parcel = 45000
dm['tc'] = transaction_cost_per_parcel
dm['tc+ce_cost_pred_usd'] = (dm['ce_cost_pred_usd']+dm['tc'])
dm['tc+ce_cost_$/tco2'] = (dm['tc+ce_cost_pred_usd'] / dm['tco2_loss_pred_100yr'])
dm['ce_cost_pred_usd_sum_with_transaction'] = (dm['tc+ce_cost_pred_usd']).cumsum()
dm_sorted5 = dm.sort_values(by='tc+ce_cost_$/tco2')
dm_sorted5


i_scenario_5 = dm_sorted5['ce_cost_pred_usd_sum_with_transaction'].lt(BUDGET)
selected_parcels_scenario_5 = dm_sorted5[i_scenario_5]
n_scenario_5 = len(selected_parcels_scenario_5)
ha_pc_scenario_5 = selected_parcels_scenario_5['ha'].sum()
ha_loss_avo_scenario_5 = selected_parcels_scenario_5['ha_f_loss_pred_100yr'].sum()
tco2_avo_scenario_5 = selected_parcels_scenario_5['tco2_loss_pred_100yr'].sum()
budget_scenario_5 = selected_parcels_scenario_5['tc+ce_cost_pred_usd'].sum()
avg_usd_tco2_scenario_5 = budget_scenario_5 / tco2_avo_scenario_5

scenario_5_values = {
    'Scenario': 'Scenario 5',
    'n': n_scenario_5,
    'ha_pc': ha_pc_scenario_5,
    'ha_loss_avo': ha_loss_avo_scenario_5,
    'tco2_avo': tco2_avo_scenario_5,
    'budget': budget_scenario_5,
    'avg_usd_tco2': avg_usd_tco2_scenario_5
}
#pd.DataFrame([scenario_5_values]).to_csv("/Users/andywang/desktop/2_carbon_prices_summary.csv", index=False, mode='a', header=False)


fig5, ax5 = plt.subplots(figsize=(17, 10))
shp_ma_outline.plot(ax=ax5, color='none', edgecolor='black')
selected_parcels_scenario_5_gdf = shp_c.join(selected_parcels_scenario_5, on='pid')
selected_parcels_scenario_5_gdf.plot(ax=ax5, markersize='ha', color='blue', alpha=0.7)
ax5.axis('off')
ax5.set_title("Scenario 5", fontsize=18)
#plt.savefig("/Users/andywang/desktop/2_selected_parcels_scenario_5.png", dpi=200)

Text(0.5, 1.0, 'Scenario 5')


random_noise = np.random.normal(scale=0.25, size=len(dm))
dm['ls_pred_wta'] = dm['ls_pred'] + random_noise + 0.1
dm['wta_ce_cost_pred_usd'] = np.exp(dm['ls_pred_wta']) * dm['ha'] * 0.4
dm_wta = dm[dm['wta_ce_cost_pred_usd'] < dm['ce_cost_pred_usd']]
dm_wta_sorted = dm_wta.sort_values(by='ce_cost_$/tco2')
dm_wta_sorted['ce_cost_pred_usd_sum'] = dm_wta_sorted['ce_cost_pred_usd'].cumsum()
i_scenario_6 = dm_wta_sorted['ce_cost_pred_usd_sum'].lt(BUDGET)


selected_parcels_scenario_6 = dm_wta_sorted[i_scenario_6]
n_scenario_6 = len(selected_parcels_scenario_6)
ha_pc_scenario_6 = selected_parcels_scenario_6['ha'].sum()
ha_loss_avo_scenario_6 = selected_parcels_scenario_6['ha_f_loss_pred_100yr'].sum()
tco2_avo_scenario_6 = selected_parcels_scenario_6['tco2_loss_pred_100yr'].sum()
budget_scenario_6 = selected_parcels_scenario_6['ce_cost_pred_usd'].sum()
avg_usd_tco2_scenario_6 = budget_scenario_6 / tco2_avo_scenario_6

scenario_6_values = {
    'Scenario': 'Scenario 6',
    'n': n_scenario_6,
    'ha_pc': ha_pc_scenario_6,
    'ha_loss_avo': ha_loss_avo_scenario_6,
    'tco2_avo': tco2_avo_scenario_6,
    'budget': budget_scenario_6,
    'avg_usd_tco2': avg_usd_tco2_scenario_6
}
#pd.DataFrame([scenario_6_values]).to_csv("/Users/andywang/desktop/2_carbon_prices_summary.csv", index=False, mode='a', header=False)


n_scenario_6 = len(selected_parcels_scenario_6)
ha_pc_scenario_6 = selected_parcels_scenario_6['ha'].sum()
ha_loss_avo_scenario_6 = selected_parcels_scenario_6['ha_f_loss_pred_100yr'].sum()
tco2_avo_scenario_6 = selected_parcels_scenario_6['tco2_loss_pred_100yr'].sum()
budget_scenario_6 = selected_parcels_scenario_6['ce_cost_pred_usd'].sum()
avg_usd_tco2_scenario_6 = budget_scenario_6 / tco2_avo_scenario_6

scenario_6_values = {
    'Scenario': 'Scenario 6',
    'n': n_scenario_6,
    'ha_pc': ha_pc_scenario_6,
    'ha_loss_avo': ha_loss_avo_scenario_6,
    'tco2_avo': tco2_avo_scenario_6,
    'budget': budget_scenario_6,
    'avg_usd_tco2': avg_usd_tco2_scenario_6
}
#pd.DataFrame([scenario_6_values]).to_csv("/Users/andywang/desktop/2_carbon_prices_summary.csv", index=False, mode='a', header=False)


fig6, ax6 = plt.subplots(figsize=(17, 10))
shp_ma_outline.plot(ax=ax6, color='none', edgecolor='black')
selected_parcels_scenario_6_gdf = shp_c.join(selected_parcels_scenario_6, on='pid')
selected_parcels_scenario_6_gdf.plot(ax=ax6, markersize='ha', color='blue', alpha=0.7)
ax6.axis('off')
ax6.set_title("Scenario 6", fontsize=18)
#plt.savefig("/Users/andywang/desktop/2_selected_parcels_scenario_6.png", dpi=200)

Text(0.5, 1.0, 'Scenario 6')


'''
Carbon offset projects are based on counterfactual thinking
'''

'\nCarbon offset projects are based on counterfactual thinking\n'


dm.head()


# Scenario 1: Optimal Allocation
dm_opt = dm.sort_values(by='ce_cost_$/tco2')
dm_opt['ce_cost_pred_usd_sum'] = dm_opt['ce_cost_pred_usd'].cumsum()
i_scenario_opt = dm_opt['ce_cost_pred_usd_sum'].lt(BUDGET*0.8)
selected_parcels_opt = dm_opt[i_scenario_opt]

# Calculating the outcomes for Scenario 1
scenario_opt_results = {
    'n': len(selected_parcels_opt),
    'ha_pc': selected_parcels_opt['ha'].sum(),
    'ha_loss_avo': selected_parcels_opt['ha_f_loss_pred_100yr'].sum(),
    'tco2_avo': selected_parcels_opt['tco2_loss_pred_100yr'].sum(),
    'budget': selected_parcels_opt['ce_cost_pred_usd'].sum(),
    'avg_usd_tco2': selected_parcels_opt['ce_cost_pred_usd'].sum() / selected_parcels_scenario_1['tco2_loss_pred_100yr'].sum()
}

# Scenario 2: Panic Buying
dm_panic = dm[dm['p_f_2010'] > 0].sort_values(by='fc_pred')
dm_panic['ce_cost_pred_usd_sum'] = dm_panic['ce_cost_pred_usd'].cumsum()
i_scenario_panic = dm_panic['ce_cost_pred_usd_sum'].lt(BUDGET*0.2)
selected_parcels_scenario_panic = dm_sorted2[i_scenario_panic]

# Calculating the outcomes for Scenario 2
scenario_panic_results = {
    'n': len(selected_parcels_scenario_panic),
    'ha_pc': selected_parcels_scenario_panic['ha'].sum(),
    'ha_loss_avo': selected_parcels_scenario_panic['ha_f_loss_pred_100yr'].sum(),
    'tco2_avo': selected_parcels_scenario_panic['tco2_loss_pred_100yr'].sum(),
    'budget': selected_parcels_scenario_panic['ce_cost_pred_usd'].sum(),
    'avg_usd_tco2': selected_parcels_scenario_panic['ce_cost_pred_usd'].sum() / selected_parcels_scenario_panic['tco2_loss_pred_100yr'].sum()
}
scenario_panic_results
#scenario_df.to_csv("/Users/andywang/desktop/2_carbon_prices_summary.csv", index=False)

#pd.DataFrame([scenario_6_values]).to_csv("/Users/andywang/desktop/2_carbon_prices_summary.csv", index=False, mode='a', header=False)

{'n': 9,
 'ha_pc': 53.827,
 'ha_loss_avo': 7.3386,
 'tco2_avo': 2935.4399999999996,
 'budget': 13500891.525943596,
 'avg_usd_tco2': 4599.273541937017}


dm['tax_offer'] = dm['acq_cost_pred_usd'] * 0.5
dm['tax_offer'] = dm['tax_offer'].clip(upper=75000)


dm['flatrate'] = dm['p_f_2010']/100 * dm['ha'] * 1210
i_scenario_tax = dm['tax_offer'] < dm['flatrate']

print(dm[i_scenario_tax]['flatrate'].sum())
selected_parcels_scenario_tax = dm[i_scenario_tax]

68928753.8709


n_scenario_tax = len(selected_parcels_scenario_tax)
ha_pc_scenario_tax = selected_parcels_scenario_tax['ha'].sum()
ha_loss_avo_scenario_tax = selected_parcels_scenario_tax['ha_f_loss_pred_100yr'].sum()
tco2_avo_scenario_tax = selected_parcels_scenario_tax['tco2_loss_pred_100yr'].sum()
budget_scenario_tax = dm[i_scenario_tax]['flatrate'].sum()
avg_usd_tco2_scenario_tax = budget_scenario_tax / tco2_avo_scenario_tax
scenario_tax_values = {
    'Scenario': 'Scenario Tax Offer',
    'n': n_scenario_tax,
    'ha_pc': ha_pc_scenario_tax,
    'ha_loss_avo': ha_loss_avo_scenario_tax,
    'tco2_avo': tco2_avo_scenario_tax,
    'budget': budget_scenario_tax,
    'avg_usd_tco2': avg_usd_tco2_scenario_tax
}
scenario_tax_values

{'Scenario': 'Scenario Tax Offer',
 'n': 587,
 'ha_pc': 66944.37700000001,
 'ha_loss_avo': 6057.395807710985,
 'tco2_avo': 2422958.323084394,
 'budget': 68928753.8709,
 'avg_usd_tco2': 28.448179737220826}


fig, ax10 = plt.subplots(figsize=(17, 10))
shp_ma_outline.plot(ax=ax10, color='none', edgecolor='black')
selected_parcels_scenario_tax_gdf = shp_c.join(selected_parcels_scenario_tax, on='pid')
selected_parcels_scenario_tax_gdf.plot(ax=ax10, markersize='ha', color='blue', alpha=0.7)
ax10.axis('off')
ax10.set_title("Tax Offer vs Florest Flat Rate", fontsize=18)

Text(0.5, 1.0, 'Tax Offer vs Florest Flat Rate')


dm['carbonrate']= dm['tco2_loss_pred_100yr'] *13.1
i_scenario_carbonrate = dm['tax_offer'] < dm['carbonrate']
print(dm[i_scenario_carbonrate]['carbonrate'].sum())

69947096.3479845


selected_parcels_scenario_carbonrate = dm[i_scenario_carbonrate]
n_scenario_carbonrate = len(selected_parcels_scenario_carbonrate)
ha_pc_scenario_carbonrate = selected_parcels_scenario_carbonrate['ha'].sum()
ha_loss_avo_scenario_carbonrate = selected_parcels_scenario_carbonrate['ha_f_loss_pred_100yr'].sum()
tco2_avo_scenario_carbonrate = selected_parcels_scenario_carbonrate['tco2_loss_pred_100yr'].sum()
budget_scenario_carbonrate = dm[i_scenario_carbonrate]['carbonrate'].sum()
avg_usd_tco2_scenario_carbonrate = budget_scenario_carbonrate / tco2_avo_scenario_carbonrate
scenario_carbonrate_values = {
    'Scenario': 'Scenario Tax Offer vs Carbon Flat Rate',
    'n': n_scenario_carbonrate,
    'ha_pc': ha_pc_scenario_carbonrate,
    'ha_loss_avo': ha_loss_avo_scenario_carbonrate,
    'tco2_avo': tco2_avo_scenario_carbonrate,
    'budget': budget_scenario_carbonrate,
    'avg_usd_tco2': avg_usd_tco2_scenario_carbonrate
}
scenario_carbonrate_values

{'Scenario': 'Scenario Tax Offer vs Carbon Flat Rate',
 'n': 555,
 'ha_pc': 43257.757,
 'ha_loss_avo': 13348.682509157348,
 'tco2_avo': 5339473.003662939,
 'budget': 69947096.3479845,
 'avg_usd_tco2': 13.099999999999998}


fig, ax11 = plt.subplots(figsize=(17, 10))
shp_ma_outline.plot(ax=ax11, color='none', edgecolor='black')
selected_parcels_scenario_carbonrate_gdf = shp_c.join(selected_parcels_scenario_carbonrate, on='pid')
selected_parcels_scenario_carbonrate_gdf.plot(ax=ax11, markersize='ha', color='blue', alpha=0.7)
ax11.axis('off')
ax11.set_title("Tax Offer vs Carbon Flat Rate", fontsize=18)

Text(0.5, 1.0, 'Tax Offer vs Carbon Flat Rate')


dm_sorted_opt = dm.sort_values(by='ce_cost_$/tco2')
dm_sorted_opt['ce_cost_pred_usd_sum'] = dm_sorted_opt['ce_cost_pred_usd'].cumsum()
dm_sorted_pan = dm[dm['p_f_2010'] > 0].sort_values(by='fc_pred')
dm_sorted_pan['ce_cost_pred_usd_sum'] = dm_sorted_pan['ce_cost_pred_usd'].cumsum()


def simulate_allocation(budget_allocation, dm_sorted_opt, dm_sorted_pan, BUDGET):
    """
    Simulates the allocation of budget between the two scenarios and calculates the total CO2 avoided 
    and average cost per tCO2 for the combined allocation.
    """
    # Allocating budget according to the given ratio
    budget_scenario_opt = BUDGET * budget_allocation
    budget_scenario_pan = BUDGET * (1 - budget_allocation)

    # Scenario 1
    i_scenario_opt = dm_sorted_opt['ce_cost_pred_usd_sum'].lt(budget_scenario_opt)
    selected_scenario_opt = dm_sorted_opt[i_scenario_opt]

    # Scenario 2
    i_scenario_pan = dm_sorted_pan['ce_cost_pred_usd_sum'].lt(budget_scenario_pan)
    selected_scenario_pan = dm_sorted_pan[i_scenario_pan]

    # Combined results
    combined_selected = pd.concat([selected_scenario_opt, selected_scenario_pan])
    total_co2_avoided = combined_selected['tco2_loss_pred_100yr'].sum()
    total_budget_used = combined_selected['ce_cost_pred_usd'].sum()
    avg_cost_per_tco2 = total_budget_used / total_co2_avoided if total_co2_avoided > 0 else np.inf

    return total_co2_avoided, avg_cost_per_tco2

allocations = np.arange(0, 1.05, 0.05)  # Allocating from 0% to 100% in 5% increments for Scenario 1
results = [simulate_allocation(allocation, dm_sorted_opt, dm_sorted_pan, BUDGET) for allocation in allocations]

# Pairing results with allocations
results_with_allocations = list(zip(allocations, results))
results_with_allocations

[(0.0, (18688.224000000002, 3661.831554417337)),
 (0.05, (208408.97933571535, 323.76859366971104)),
 (0.1, (349440.7591711177, 199.9041432084736)),
 (0.15000000000000002, (472187.6891205429, 139.24262811426237)),
 (0.2, (590313.2495407679, 117.41693307236335)),
 (0.25, (699898.3108133498, 93.16154606766555)),
 (0.30000000000000004, (803269.6498168104, 85.4234011409276)),
 (0.35000000000000003, (904628.2792642632, 52.2663046392076)),
 (0.4, (1001494.9751070074, 50.588448312530694)),
 (0.45, (1105043.2615579497, 49.196302369869656)),
 (0.5, (1198860.8999801083, 48.22000418811502)),
 (0.55, (1293750.6744031503, 47.432664764532504)),
 (0.6000000000000001, (1380677.2648041507, 46.857144717419835)),
 (0.65, (1472461.0603783657, 46.371398859145835)),
 (0.7000000000000001, (1518150.9516055286, 44.7606240096028)),
 (0.75, (1641723.4929221836, 42.00639230733035)),
 (0.8, (1728889.8029860104, 40.19401911022835)),
 (0.8500000000000001, (1811077.036106914, 37.253067114854645)),
 (0.9, (1893917.2093978468, 35.08576656362387)),
 (0.9500000000000001, (1974863.0185223534, 35.40127655152479)),
 (1.0, (2055343.026564267, 34.05601719336201))]


# Extracting total CO2 avoided and average cost per tCO2 for each allocation
total_co2_avoided = [result[0] for _, result in results_with_allocations]
avg_cost_per_tco2 = [result[1] for _, result in results_with_allocations]

# Plotting the results
fig, ax1 = plt.subplots()

# Bar graph for total CO2 avoided
ax1.bar(allocations, total_co2_avoided, width=0.03, color='b', align='center')
ax1.set_xlabel('Allocation to Scenario 1')
ax1.set_ylabel('Total CO2 Avoided (tCO2)', color='b')
ax1.tick_params('y', colors='b')

# Creating a second y-axis for average cost per tCO2
ax2 = ax1.twinx()
ax2.plot(allocations, avg_cost_per_tco2, 'r-')
ax2.set_ylabel('Average Cost per tCO2 (USD)', color='r')
ax2.tick_params('y', colors='r')

plt.title('Total CO2 Avoided and Average Cost per tCO2 by Allocation')
plt.show()


def simulate_allocation2(budget_allocation, dm_sorted_opt, dm_sorted_pan, BUDGET):
    """
    Simulates the allocation of budget between the two scenarios and calculates the total CO2 avoided 
    and average cost per tCO2 for the combined allocation.
    """
    # Allocating budget according to the given ratio
    budget_scenario_opt = BUDGET * budget_allocation
    budget_scenario_pan = BUDGET * (1 - budget_allocation)

    # Scenario 1
    i_scenario_opt = dm_sorted_opt['ce_cost_pred_usd_sum'].lt(budget_scenario_opt)
    selected_scenario_opt = dm_sorted_opt[i_scenario_opt]

    # Scenario 2
    i_scenario_pan = dm_sorted_pan['ce_cost_pred_usd_sum'].lt(budget_scenario_pan)
    selected_scenario_pan = dm_sorted_pan[i_scenario_pan]

    # Combined results
    combined_selected = pd.concat([selected_scenario_opt, selected_scenario_pan])
    total_co2_avoided = (selected_scenario_opt['tco2_loss_pred_100yr'].sum()) + ((selected_scenario_pan['tco2_loss_pred_100yr']*15).sum())
    total_budget_used = combined_selected['ce_cost_pred_usd'].sum()
    avg_cost_per_tco2 = total_budget_used / total_co2_avoided if total_co2_avoided > 0 else np.inf

    return total_co2_avoided, avg_cost_per_tco2

allocations = np.arange(0, 1.05, 0.001)  # Allocating from 0% to 100% in 5% increments for Scenario 1
results2 = [simulate_allocation2(allocation, dm_sorted_opt, dm_sorted_pan, BUDGET) for allocation in allocations]

# Pairing results with allocations
results_with_allocations2 = list(zip(allocations, results2))
results_with_allocations2

[(0.0, (280323.36000000004, 244.12210362782247)),
 (0.001, (288713.8439679416, 237.26656636083797)),
 (0.002, (294996.2639833978, 232.4384172849466)),
 (0.003, (300596.9815912618, 228.32142810910585)),
 (0.004, (300596.9815912618, 228.32142810910585)),
 (0.005, (309825.13321610674, 221.88929843743156)),
 (0.006, (316876.0957418652, 217.24652185580013)),
 (0.007, (316876.0957418652, 217.24652185580013)),
 (0.008, (327729.0138672049, 210.51490343465036)),
 (0.009000000000000001, (331758.5987963735, 208.13088225289147)),
 (0.01, (331758.5987963735, 208.13088225289147)),
 (0.011, (341425.6867963736, 202.6468853055562)),
 (0.012, (344588.328306393, 200.9237092670085)),
 (0.013000000000000001, (349668.3489565756, 198.22663018346896)),
 (0.014, (355837.8782911397, 195.06133766776333)),
 (0.015, (359600.0045853751, 193.1934664265863)),
 (0.016, (363904.37898900127, 191.109884838842)),
 (0.017, (367960.4945795689, 189.1992477427866)),
 (0.018000000000000002, (369517.31664440315, 188.47967546312887)),
 (0.019, (373873.61365684855, 186.50146869741778)),
 (0.02, (378433.83365830104, 184.48199597444028)),
 (0.021, (382624.0307919561, 182.67188575369082)),
 (0.022, (386095.6373751653, 181.2055838128344)),
 (0.023, (383237.2373751653, 182.13912526566745)),
 (0.024, (389678.0340913662, 179.45627608916558)),
 (0.025, (392604.841010284, 177.50995823357292)),
 (0.026000000000000002, (396479.51006050914, 175.97124456981916)),
 (0.027, (400538.6202903059, 174.39338756510426)),
 (0.028, (403779.022465399, 173.15779819657237)),
 (0.029, (407463.7217820729, 171.77897739117378)),
 (0.03, (405849.34711978433, 172.0418983990292)),
 (0.031, (409691.6611566987, 170.62518099579214)),
 (0.032, (352433.77304579003, 197.11864800168212)),
 (0.033, (355572.4123710933, 195.56562187792233)),
 (0.034, (355572.4123710933, 195.56562187792233)),
 (0.035, (364465.4886499952, 191.31453042635528)),
 (0.036000000000000004, (366603.39786894445, 190.3250935284417)),
 (0.037, (370574.2855161191, 188.51859195385242)),
 (0.038, (372572.11358663917, 187.62462469407663)),
 (0.039, (373252.25384664215, 178.89315529807138)),
 (0.04, (373252.25384664215, 178.89315529807138)),
 (0.041, (379784.83542128076, 176.1953143073245)),
 (0.042, (380590.11700571456, 175.86930235084364)),
 (0.043000000000000003, (384844.6121527576, 174.1705963298739)),
 (0.044, (387339.25902592135, 173.1924496877539)),
 (0.045, (392393.63683757844, 171.2493666456154)),
 (0.046, (394744.16211920226, 170.3628569292187)),
 (0.047, (397309.9401704534, 169.40770484684194)),
 (0.048, (401084.66194747627, 168.0277480724304)),
 (0.049, (404695.9233357154, 166.73328851817837)),
 (0.05, (404695.9233357154, 166.73328851817837)),
 (0.051000000000000004, (410752.6815293126, 164.61526396354935)),
 (0.052000000000000005, (413017.0259919108, 163.83992612951985)),
 (0.053, (415284.6237039282, 163.07261079617254)),
 (0.054, (419802.61862060556, 161.56956205651335)),
 (0.055, (422404.58482978435, 160.7203730202542)),
 (0.056, (425374.28873072285, 159.76441231572184)),
 (0.057, (428125.21038101305, 158.89117334349015)),
 (0.058, (429877.95166145643, 158.34139604721594)),
 (0.059000000000000004, (433226.24822982715, 157.3037619702776)),
 (0.06, (436960.37935113243, 156.1664978479264)),
 (0.061, (439018.2236581172, 155.54862184847946)),
 (0.062, (439018.2236581172, 155.54862184847946)),
 (0.063, (439018.2236581172, 155.54862184847946)),
 (0.064, (447566.51517171843, 153.04481546025642)),
 (0.065, (450845.5290900913, 152.1100597542526)),
 (0.066, (453454.5242858703, 151.3762939704835)),
 (0.067, (453454.5242858703, 151.3762939704835)),
 (0.068, (457892.50090475427, 150.14751080671147)),
 (0.069, (461388.72695372836, 149.19665907402444)),
 (0.07, (465719.8983784538, 148.03872388287684)),
 (0.07100000000000001, (467501.8613122781, 147.56881998365057)),
 (0.07200000000000001, (470792.55895586975, 146.7105392463034)),
 (0.073, (472007.24949469394, 146.39681257402617)),
 (0.074, (476291.7734618398, 145.303819722673)),
 (0.075, (478230.08559919096, 144.8157929962546)),
 (0.076, (482567.0230449291, 143.73881748944953)),
 (0.077, (482567.0230449291, 143.73881748944953)),
 (0.078, (486488.62774646626, 142.78315777850707)),
 (0.079, (490318.4563113604, 141.86546920231984)),
 (0.08, (493528.78760656324, 141.10748297354905)),
 (0.081, (495739.19889395055, 140.5914506163418)),
 (0.082, (495739.19889395055, 140.5914506163418)),
 (0.083, (495739.19889395055, 140.5914506163418)),
 (0.084, (504869.6820776064, 138.50791285229096)),
 (0.085, (506316.8146928262, 138.1848779785772)),
 (0.08600000000000001, (500013.96579206135, 138.40637202091614)),
 (0.08700000000000001, (502363.68156453595, 137.87884078411696)),
 (0.088, (504409.4666319285, 137.42394059212313)),
 (0.089, (507703.79646929255, 136.69984713311092)),
 (0.09, (511160.47404613445, 135.9507364566432)),
 (0.091, (511160.47404613445, 135.9507364566432)),
 (0.092, (511160.47404613445, 135.9507364566432)),
 (0.093, (519418.33873952867, 134.20253908016448)),
 (0.094, (521588.3468784452, 133.7525741856015)),
 (0.095, (523994.61129577603, 133.2581373110374)),
 (0.096, (525423.6396596634, 132.96666880833)),
 (0.097, (529993.002688166, 132.04589694736475)),
 (0.098, (528735.7809496236, 131.91767334868416)),
 (0.099, (531113.5265943466, 131.4452014856597)),
 (0.1, (532717.5591711176, 131.12887751046026)),
 (0.101, (534843.8527912907, 130.71250884003157)),
 (0.10200000000000001, (536703.5473733388, 129.57211955185755)),
 (0.10300000000000001, (536890.9974453915, 129.536136353385)),
 (0.10400000000000001, (541330.3946336457, 128.69150529805495)),
 (0.105, (541330.3946336457, 128.69150529805495)),
 (0.106, (545760.5855520514, 127.86255053857808)),
 (0.107, (549874.4043234477, 127.10525227415737)),
 (0.108, (552234.5218079343, 126.67647113548112)),
 (0.109, (533160.5680381829, 118.0557292627597)),
 (0.11, (534164.3349295823, 117.88497526442607)),
 (0.111, (537374.247575688, 117.34359258561803)),
 (0.112, (539462.0340519447, 116.99523335522038)),
 (0.113, (539462.0340519447, 116.99523335522038)),
 (0.114, (544480.3200676299, 116.1690487730297)),
 (0.115, (544480.3200676299, 116.1690487730297)),
 (0.116, (544480.3200676299, 116.1690487730297)),
 (0.117, (553638.8633606753, 114.70097611530228)),
 (0.11800000000000001, (553638.8633606753, 114.70097611530228)),
 (0.11900000000000001, (558650.8394436841, 113.91852875624943)),
 (0.12, (558650.8394436841, 113.91852875624943)),
 (0.121, (563834.0327637601, 113.12477221197025)),
 (0.122, (565734.6757205155, 112.83751054717287)),
 (0.123, (568501.0164592904, 112.42319270919165)),
 (0.124, (571343.3848273718, 112.00216559106583)),
 (0.125, (573183.017854233, 111.73205273722888)),
 (0.126, (576490.5459604702, 111.25101503106767)),
 (0.127, (577650.4205263024, 111.08372880180481)),
 (0.128, (581043.2146936544, 110.59859863280779)),
 (0.129, (583318.3695082249, 110.27670640647764)),
 (0.13, (586503.939959044, 109.83104677365)),
 (0.131, (587839.7878977822, 109.64571942507959)),
 (0.132, (587839.7878977822, 109.64571942507959)),
 (0.133, (593874.8641012921, 108.81969396651516)),
 (0.134, (596062.7013105411, 108.52476076996172)),
 (0.135, (598803.8537820401, 108.15869477381833)),
 (0.136, (599097.7907706681, 108.11973601490624)),
 (0.137, (599097.7907706681, 108.11973601490624)),
 (0.138, (599097.7907706681, 108.11973601490624)),
 (0.139, (599097.7907706681, 108.11973601490624)),
 (0.14, (610976.0762745829, 106.57787536706317)),
 (0.14100000000000001, (613282.9938773179, 106.2856983026789)),
 (0.14200000000000002, (614929.1643524455, 106.0786268656114)),
 (0.14300000000000002, (616977.2231716219, 105.82266920156547)),
 (0.14400000000000002, (619693.7303328018, 105.4861570104073)),
 (0.145, (619693.7303328018, 105.4861570104073)),
 (0.146, (624825.6268319181, 104.85868096236175)),
 (0.147, (624825.6268319181, 104.85868096236175)),
 (0.148, (628378.3644604621, 104.43031418155745)),
 (0.149, (632306.9070855033, 103.96241960929284)),
 (0.15, (632735.7691205429, 103.91170849678781)),
 (0.151, (636984.6963416561, 103.41345287018373)),
 (0.152, (639642.8853156375, 103.10525448108808)),
 (0.153, (641992.3274013412, 102.83507379228585)),
 (0.154, (642987.1205005466, 102.72134616239529)),
 (0.155, (646425.7080968248, 102.3314365361484)),
 (0.156, (649398.3208933763, 101.99809797815621)),
 (0.157, (649398.3208933763, 101.99809797815621)),
 (0.158, (653411.1464325943, 101.55317085224753)),
 (0.159, (656330.2862286085, 101.23347619912903)),
 (0.16, (658471.8407853036, 101.00075791454245)),
 (0.161, (660293.4520971515, 100.80412016211739)),
 (0.162, (660293.4520971515, 100.80412016211739)),
 (0.163, (664907.5606641974, 100.31100511905426)),
 (0.164, (668105.1827653403, 99.9733533101342)),
 (0.165, (669928.0652482576, 99.7824125322852)),
 (0.166, (671960.5048018614, 99.57093906815477)),
 (0.167, (674743.9034732745, 99.28369596995266)),
 (0.168, (677395.2729331465, 99.01247044877637)),
 (0.169, (679148.7653936103, 98.8344602705621)),
 (0.17, (682233.0495026774, 98.52380489993566)),
 (0.171, (684413.103045214, 98.30603265122684)),
 (0.17200000000000001, (685379.6518221999, 98.20993049014673)),
 (0.17300000000000001, (688935.1996898229, 97.85886046557499)),
 (0.17400000000000002, (690114.8340929761, 97.74323664249005)),
 (0.17500000000000002, (693083.8552168652, 97.45428383643804)),
 (0.176, (694174.6179645985, 97.34884802025516)),
 (0.177, (694174.6179645985, 97.34884802025516)),
 (0.178, (700369.5121673498, 96.75633131508656)),
 (0.179, (700369.5121673498, 96.75633131508656)),
 (0.18, (703362.4245648722, 96.47403644145813)),
 (0.181, (707144.26660278, 96.12091256249995)),
 (0.182, (709544.2544063726, 95.8989019879502)),
 (0.183, (709544.2544063726, 95.8989019879502)),
 (0.184, (712784.2661111713, 95.60168000061336)),
 (0.185, (712784.2661111713, 95.60168000061336)),
 (0.186, (717875.7683000085, 95.14031459056804)),
 (0.187, (717875.7683000085, 95.14031459056804)),
 (0.188, (717875.7683000085, 95.14031459056804)),
 (0.189, (725056.8748309747, 94.50101248278338)),
 (0.19, (727526.1589984655, 94.28413181243603)),
 (0.191, (729888.3709504185, 94.07835601052633)),
 (0.192, (729888.3709504185, 94.07835601052633)),
 (0.193, (733429.4541148982, 93.77243283171589)),
 (0.194, (736808.9638116169, 93.48332372352617)),
 (0.195, (736808.9638116169, 93.48332372352617)),
 (0.196, (740890.4261470108, 93.13785354576332)),
 (0.197, (742461.3342627496, 93.00594400151391)),
 (0.198, (744529.5018685127, 92.83326807717539)),
 (0.199, (748319.6513457019, 92.51945230889893)),
 (0.2, (750861.3295407679, 92.31101481208233)),
 (0.201, (751999.2416194648, 92.218171297144)),
 (0.202, (754128.0571376926, 92.04527101335131)),
 (0.203, (754128.0571376926, 92.04527101335131)),
 (0.20400000000000001, (758340.8973684066, 91.70610990111467)),
 (0.20500000000000002, (761358.6716209534, 91.46557866707401)),
 (0.20600000000000002, (763355.1834041034, 91.30758103564335)),
 (0.20700000000000002, (766297.167251514, 91.0764761415436)),
 (0.20800000000000002, (766297.167251514, 91.0764761415436)),
 (0.209, (770134.1119280759, 90.77834257497054)),
 (0.21, (769913.3172584383, 90.68768548036982)),
 (0.211, (772201.165660851, 90.51185702382942)),
 (0.212, (774422.7223312644, 90.34229146227237)),
 (0.213, (775325.2001757037, 89.42555530237674)),
 (0.214, (775325.2001757037, 89.42555530237674)),
 (0.215, (775325.2001757037, 89.42555530237674)),
 (0.216, (781235.1653724734, 88.98760113797827)),
 (0.217, (781235.1653724734, 88.98760113797827)),
 (0.218, (781235.1653724734, 88.98760113797827)),
 (0.219, (788354.7364356115, 88.46878250107159)),
 (0.22, (790945.2202629738, 88.28236301729834)),
 (0.221, (793406.1285671224, 88.10648428416984)),
 (0.222, (794739.7702038572, 88.01167560135569)),
 (0.223, (792643.1302038572, 88.01226394510375)),
 (0.224, (797619.2540121548, 87.66054172926613)),
 (0.225, (795255.8800965825, 79.71911136987276)),
 (0.226, (798326.5675279559, 79.53433494785655)),
 (0.227, (801566.5563707317, 79.34109484158655)),
 (0.228, (802244.893170045, 79.30084546764832)),
 (0.229, (804345.2086027751, 79.1767193516618)),
 (0.23, (806955.3557080526, 79.02342336220049)),
 (0.231, (806955.3557080526, 79.02342336220049)),
 (0.232, (812468.4081362969, 78.7033240794976)),
 (0.233, (813858.901621398, 78.62337373871996)),
 (0.234, (817042.1817340539, 78.44154200969506)),
 (0.23500000000000001, (818959.9726436192, 78.33280682401968)),
 (0.23600000000000002, (821275.2867328429, 78.20226493515555)),
 (0.23700000000000002, (822139.6466072472, 78.15374256485434)),
 (0.23800000000000002, (825666.2131558554, 77.95696918680822)),
 (0.23900000000000002, (827332.8946999549, 77.86469687014052)),
 (0.24, (828729.0380940626, 77.78774734111815)),
 (0.241, (831866.8581510512, 77.6159311054035)),
 (0.242, (834083.4293254395, 77.49556047366403)),
 (0.243, (835795.9286086089, 77.40314622842011)),
 (0.244, (837896.0994257838, 77.29040107324549)),
 (0.245, (840472.9301223276, 77.15295835939146)),
 (0.246, (842824.2854294993, 77.02828130442657)),
 (0.247, (842824.2854294993, 77.02828130442657)),
 (0.248, (842824.2854294993, 77.02828130442657)),
 (0.249, (842824.2854294993, 77.02828130442657)),
 (0.25, (851519.0388133498, 76.57328345398466)),
 (0.251, (852783.609417615, 76.50791799487546)),
 (0.252, (854863.867832747, 76.40084492998572)),
 (0.253, (856905.1848165947, 76.29633738439829)),
 (0.254, (856905.1848165947, 76.29633738439829)),
 (0.255, (856905.1848165947, 76.29633738439829)),
 (0.256, (863997.0536635033, 75.93723967245062)),
 (0.257, (863997.0536635033, 75.93723967245062)),
 (0.258, (863997.0536635033, 75.93723967245062)),
 (0.259, (863997.0536635033, 75.93723967245062)),
 (0.26, (871899.8238389111, 75.54452699187877)),
 (0.261, (871899.8238389111, 75.54452699187877)),
 (0.262, (871899.8238389111, 75.54452699187877)),
 (0.263, (878642.0858993162, 75.2151126125004)),
 (0.264, (881578.7291506311, 75.07380129472081)),
 (0.265, (881578.7291506311, 75.07380129472081)),
 (0.266, (881578.7291506311, 75.07380129472081)),
 (0.267, (887468.9493053427, 74.79350247367438)),
 (0.268, (889640.2700941954, 74.69130284775541)),
 (0.269, (892358.3720924594, 74.56435966434904)),
 (0.27, (892358.3720924594, 74.56435966434904)),
 (0.271, (896472.5219775252, 74.37380704562877)),
 (0.272, (897066.0552528873, 74.34647085563816)),
 (0.273, (899306.3332382625, 74.24363565712781)),
 (0.274, (899306.3332382625, 74.24363565712781)),
 (0.275, (899306.3332382625, 74.24363565712781)),
 (0.276, (906987.1030810544, 73.89531502665828)),
 (0.277, (907736.9387018699, 73.86165044762157)),
 (0.278, (907736.9387018699, 73.86165044762157)),
 (0.279, (907736.9387018699, 73.86165044762157)),
 (0.28, (907736.9387018699, 73.86165044762157)),
 (0.281, (917072.9896207985, 73.44732587841968)),
 (0.28200000000000003, (917072.9896207985, 73.44732587841968)),
 (0.28300000000000003, (917072.9896207985, 73.44732587841968)),
 (0.28400000000000003, (924023.064850831, 73.1444631343147)),
 (0.28500000000000003, (924349.6675257566, 73.1303542960508)),
 (0.28600000000000003, (926685.9053248056, 73.02982218715236)),
 (0.28700000000000003, (926685.9053248056, 73.02982218715236)),
 (0.28800000000000003, (932145.8147621285, 72.79708381582655)),
 (0.289, (933578.0553352935, 72.73649013238114)),
 (0.29, (936432.8150581808, 72.61635569334952)),
 (0.291, (938031.0131453939, 72.54947245610957)),
 (0.292, (938031.0131453939, 72.54947245610957)),
 (0.293, (942358.8107445145, 72.3696510354346)),
 (0.294, (944740.823092622, 72.27147678960762)),
 (0.295, (946851.9346267383, 72.1848903530583)),
 (0.296, (948409.5470732339, 72.12132562293677)),
 (0.297, (950659.092952385, 72.03000102693508)),
 (0.298, (951421.1694201616, 71.99918594030302)),
 (0.299, (954890.3778168105, 71.85958421480476)),
 (0.3, (954890.3778168105, 71.85958421480476)),
 (0.301, (958267.082524488, 71.72478344675453)),
 (0.302, (958267.082524488, 71.72478344675453)),
 (0.303, (962160.1305195732, 71.5705451718161)),
 (0.304, (962160.1305195732, 71.5705451718161)),
 (0.305, (962160.1305195732, 71.5705451718161)),
 (0.306, (962160.1305195732, 71.5705451718161)),
 (0.307, (970612.0465799898, 71.24038459230063)),
 (0.308, (973026.1265799899, 71.14715051425618)),
 (0.309, (975159.1229433364, 71.06516227593252)),
 (0.31, (976631.6964190034, 71.00880351556016)),
 (0.311, (980371.0870102594, 70.86650968365086)),
 (0.312, (980921.1244758915, 70.8456955844089)),
 (0.313, (983028.6909914156, 70.7661747394687)),
 (0.314, (986462.3183934498, 70.63739275235534)),
 (0.315, (986462.3183934498, 70.63739275235534)),
 (0.316, (990675.6038521358, 70.48073737547618)),
 (0.317, (992646.962734439, 70.40799103223117)),
 (0.318, (992646.962734439, 70.40799103223117)),
 (0.319, (979591.2951021367, 46.14466033721565)),
 (0.32, (979591.2951021367, 46.14466033721565)),
 (0.321, (983079.6384156563, 46.10161224227952)),
 (0.322, (983079.6384156563, 46.10161224227952)),
 (0.323, (986698.5768213919, 46.05728731830655)),
 (0.324, (988446.9855275659, 46.03602093971247)),
 (0.325, (990899.7425428599, 46.00633991115474)),
 (0.326, (993607.0871734044, 45.97378401361214)),
 (0.327, (993607.0871734044, 45.97378401361214)),
 (0.328, (993607.0871734044, 45.97378401361214)),
 (0.329, (993607.0871734044, 45.97378401361214)),
 (0.33, (993607.0871734044, 45.97378401361214)),
 (0.331, (993607.0871734044, 45.97378401361214)),
 (0.332, (993607.0871734044, 45.97378401361214)),
 (0.333, (1007870.6387927651, 45.807357192879465)),
 (0.334, (1009534.417752848, 45.78828022713917)),
 (0.335, (1009534.417752848, 45.78828022713917)),
 (0.336, (1009534.417752848, 45.78828022713917)),
 (0.337, (1016643.5827728522, 45.70761672267661)),
 (0.338, (1018628.602212385, 45.68544575832863)),
 (0.339, (1020849.3079159656, 45.660835648365115)),
 (0.34, (1020849.3079159656, 45.660835648365115)),
 (0.341, (1023340.4970658841, 45.63339816297698)),
 (0.342, (1025100.5569352238, 45.6141198520125)),
 (0.343, (1025100.5569352238, 45.6141198520125)),
 (0.34400000000000003, (1025100.5569352238, 45.6141198520125)),
 (0.34500000000000003, (1025100.5569352238, 45.6141198520125)),
 (0.34600000000000003, (1034791.6597035888, 45.50932037801288)),
 (0.34700000000000003, (1036008.8061160499, 45.49633947704848)),
 (0.34800000000000003, (1037317.166614573, 45.48250967065748)),
 (0.34900000000000003, (1040263.1352642633, 45.451555117597216)),
 (0.35000000000000003, (1040263.1352642633, 45.451555117597216)),
 (0.35100000000000003, (1040263.1352642633, 45.451555117597216)),
 (0.352, (1040263.1352642633, 45.451555117597216)),
 (0.353, (1048922.6817314504, 45.36168629993604)),
 (0.354, (1050667.817586246, 45.34390434602361)),
 (0.355, (1051763.152508797, 45.33279027867306)),
 (0.356, (1051763.152508797, 45.33279027867306)),
 (0.357, (1051763.152508797, 45.33279027867306)),
 (0.358, (1058104.3306244076, 45.26907841688496)),
 (0.359, (1058104.3306244076, 45.26907841688496)),
 (0.36, (1058104.3306244076, 45.26907841688496)),
 (0.361, (1058104.3306244076, 45.26907841688496)),
 (0.362, (1066227.379428694, 45.18898383462343)),
 (0.363, (1068674.9943617187, 45.16514680669464)),
 (0.364, (1068674.9943617187, 45.16514680669464)),
 (0.365, (1068674.9943617187, 45.16514680669464)),
 (0.366, (1075029.8515041443, 45.103932163778545)),
 (0.367, (1075029.8515041443, 45.103932163778545)),
 (0.368, (1078641.8515041443, 45.069552403520085)),
 (0.369, (1081482.8712782052, 45.042704428904585)),
 (0.37, (1081957.1120458066, 45.03825061401699)),
 (0.371, (1081957.1120458066, 45.03825061401699)),
 (0.372, (1081957.1120458066, 45.03825061401699)),
 (0.373, (1081957.1120458066, 45.03825061401699)),
 (0.374, (1090805.5120458067, 44.956015436480875)),
 (0.375, (1090805.5120458067, 44.956015436480875)),
 (0.376, (1094702.2963103165, 44.92028289544919)),
 (0.377, (1097444.6803458408, 44.89545075971302)),
 (0.378, (1097444.6803458408, 44.89545075971302)),
 (0.379, (1097444.6803458408, 44.89545075971302)),
 (0.38, (1097444.6803458408, 44.89545075971302)),
 (0.381, (1105395.6506130889, 44.824256449122906)),
 (0.382, (1106377.3452057398, 44.815549304503634)),
 (0.383, (1108352.6735358622, 44.79807916314286)),
 (0.384, (1109708.4871474945, 44.78612678709219)),
 (0.385, (1112046.3821364664, 44.765607926589944)),
 (0.386, (1115301.3284623947, 44.73725893842764)),
 (0.387, (1117401.7742044409, 44.71912362342939)),
 (0.388, (1117401.7742044409, 44.71912362342939)),
 (0.389, (1117401.7742044409, 44.71912362342939)),
 (0.39, (1123502.9431471962, 44.66700176528907)),
 (0.391, (1125408.454372053, 44.650876656328414)),
 (0.392, (1126343.550455512, 44.64303380092687)),
 (0.393, (1129433.545948443, 44.61730196562377)),
 (0.394, (1130023.709054299, 44.61243391877287)),
 (0.395, (1132368.7603157607, 44.59316558803536)),
 (0.396, (1133936.20494226, 44.58033596572943)),
 (0.397, (1137129.8311070076, 44.55434673993035)),
 (0.398, (1137129.8311070076, 44.55434673993035)),
 (0.399, (1137129.8311070076, 44.55434673993035)),
 (0.4, (1137129.8311070076, 44.55434673993035)),
 (0.401, (1137129.8311070076, 44.55434673993035)),
 (0.402, (1137129.8311070076, 44.55434673993035)),
 (0.403, (1148952.740620424, 44.459855832893325)),
 (0.404, (1148952.740620424, 44.459855832893325)),
 (0.405, (1152519.6978856905, 44.431740909808866)),
 (0.406, (1154766.542387676, 44.414235755432216)),
 (0.40700000000000003, (1154766.542387676, 44.414235755432216)),
 (0.40800000000000003, (1159123.2668122984, 44.38050932992492)),
 (0.40900000000000003, (1160030.2438560913, 44.373541234952796)),
 (0.41000000000000003, (1160030.2438560913, 44.373541234952796)),
 (0.41100000000000003, (1164241.8147716764, 44.34132972969252)),
 (0.41200000000000003, (1167060.126905465, 44.31999478614795)),
 (0.41300000000000003, (1167766.9162008637, 44.314705772889454)),
 (0.41400000000000003, (1170305.3071389853, 44.295862040254704)),
 (0.41500000000000004, (1171485.0192665055, 44.28714182394274)),
 (0.41600000000000004, (1171485.0192665055, 44.28714182394274)),
 (0.417, (1171485.0192665055, 44.28714182394274)),
 (0.418, (1171485.0192665055, 44.28714182394274)),
 (0.419, (1180727.8294586996, 44.21949054250587)),
 (0.42, (1181179.7544269322, 44.216219737270734)),
 (0.421, (1184174.3461762918, 44.19463699689569)),
 (0.422, (1184174.3461762918, 44.19463699689569)),
 (0.423, (1184174.3461762918, 44.19463699689569)),
 (0.424, (1184174.3461762918, 44.19463699689569)),
 (0.425, (1184174.3461762918, 44.19463699689569)),
 (0.426, (1184174.3461762918, 44.19463699689569)),
 (0.427, (1184174.3461762918, 44.19463699689569)),
 (0.428, (1197475.9583915316, 44.10020823352574)),
 (0.429, (1199978.2431162565, 44.082839449599085)),
 (0.43, (1201860.9138055767, 44.06991114329871)),
 (0.431, (1204107.9480728044, 44.05463689693561)),
 (0.432, (1204640.7040083574, 44.05104580691101)),
 (0.433, (1204640.7040083574, 44.05104580691101)),
 (0.434, (1208984.4167873391, 44.02188789760931)),
 (0.435, (1210757.3620163281, 44.010050049183185)),
 (0.436, (1212837.394016328, 43.99621051371344)),
 (0.437, (1212837.394016328, 43.99621051371344)),
 (0.438, (1212837.394016328, 43.99621051371344)),
 (0.439, (1219184.3146512485, 43.95429088963534)),
 (0.44, (1220854.1659796364, 43.943340724923246)),
 (0.441, (1223334.2889821571, 43.927136134300895)),
 (0.442, (1225346.7760332383, 43.914051895862485)),
 (0.443, (1227629.3001118815, 43.89932285098565)),
 (0.444, (1229092.1788793693, 43.89002324518806)),
 (0.445, (1231498.6897075009, 43.87482387970246)),
 (0.446, (1231498.6897075009, 43.87482387970246)),
 (0.447, (1234850.7228359035, 43.85395961101607)),
 (0.448, (1236819.3197508873, 43.84176778012131)),
 (0.449, (1236819.3197508873, 43.84176778012131)),
 (0.45, (1240678.1175579496, 43.81800698991745)),
 (0.451, (1243148.3521327027, 43.80297590642868)),
 (0.452, (1244569.6761026108, 43.79443376895202)),
 (0.453, (1245252.162696369, 43.790364870139356)),
 (0.454, (1248020.884684029, 43.77391862905835)),
 (0.455, (1248020.884684029, 43.77391862905835)),
 (0.456, (1248020.884684029, 43.77391862905835)),
 (0.457, (1248020.884684029, 43.77391862905835)),
 (0.458, (1248020.884684029, 43.77391862905835)),
 (0.459, (1248020.884684029, 43.77391862905835)),
 (0.46, (1248020.884684029, 43.77391862905835)),
 (0.461, (1262132.262879228, 43.69125862023591)),
 (0.462, (1263396.0615856436, 43.683979878471106)),
 (0.463, (1265431.552214655, 43.67233202459821)),
 (0.464, (1268062.067006189, 43.657384889486856)),
 (0.465, (1268062.067006189, 43.657384889486856)),
 (0.466, (1271401.705185094, 43.6385772395376)),
 (0.467, (1271401.705185094, 43.6385772395376)),
 (0.468, (1274936.1778122378, 43.6188939478953)),
 (0.46900000000000003, (1276501.6635176167, 43.610250956735435)),
 (0.47000000000000003, (1279542.3697777884, 43.59353253522435)),
 (0.47100000000000003, (1280697.379403688, 43.58725132986373)),
 (0.47200000000000003, (1282558.0574897, 43.57718074637724)),
 (0.47300000000000003, (1284245.3855149182, 43.56811820370469)),
 (0.47400000000000003, (1286681.9183029726, 43.55518922816362)),
 (0.47500000000000003, (1289201.8775555291, 43.541917057706)),
 (0.47600000000000003, (1290844.6582439695, 43.533324846756855)),
 (0.47700000000000004, (1292282.6280260945, 43.525841656333846)),
 (0.47800000000000004, (1292282.6280260945, 43.525841656333846)),
 (0.47900000000000004, (1292282.6280260945, 43.525841656333846)),
 (0.48, (1292282.6280260945, 43.525841656333846)),
 (0.481, (1292282.6280260945, 43.525841656333846)),
 (0.482, (1301860.3077055295, 43.47676185873371)),
 (0.483, (1303231.8935995481, 43.469813812312324)),
 (0.484, (1303231.8935995481, 43.469813812312324)),
 (0.485, (1307926.5042226862, 43.44616706746732)),
 (0.486, (1309741.621709506, 43.43707619854991)),
 (0.487, (1311240.0797741916, 43.42966105106292)),
 (0.488, (1312822.0875215845, 43.42185276496366)),
 (0.489, (1314905.2476985177, 43.41162790174518)),
 (0.49, (1316711.9646005728, 43.4028151871498)),
 (0.491, (1318814.693597164, 43.392652833027235)),
 (0.492, (1318814.693597164, 43.392652833027235)),
 (0.493, (1322910.4913792717, 43.37316158440394)),
 (0.494, (1324398.507392078, 43.366133160480665)),
 (0.495, (1324398.507392078, 43.366133160480665)),
 (0.496, (1324398.507392078, 43.366133160480665)),
 (0.497, (1324398.507392078, 43.366133160480665)),
 (0.498, (1331337.7969407768, 43.333659508191154)),
 (0.499, (1334495.7559801082, 43.31904193696803)),
 (0.5, (1334495.7559801082, 43.31904193696803)),
 (0.501, (1334495.7559801082, 43.31904193696803)),
 (0.502, (1334495.7559801082, 43.31904193696803)),
 (0.503, (1334495.7559801082, 43.31904193696803)),
 (0.504, (1334495.7559801082, 43.31904193696803)),
 (0.505, (1334495.7559801082, 43.31904193696803)),
 (0.506, (1347502.5416236625, 43.259704721766575)),
 (0.507, (1348141.2465956733, 43.25683613058865)),
 (0.508, (1350630.4162415706, 43.24574282233502)),
 (0.509, (1350630.4162415706, 43.24574282233502)),
 (0.51, (1350630.4162415706, 43.24574282233502)),
 (0.511, (1356774.4678275536, 43.21859676444767)),
 (0.512, (1358932.271675483, 43.20917065001594)),
 (0.513, (1360181.836540783, 43.20372703317565)),
 (0.514, (1360181.836540783, 43.20372703317565)),
 (0.515, (1364580.1776882664, 43.184666089156984)),
 (0.516, (1364580.1776882664, 43.184666089156984)),
 (0.517, (1367934.678682162, 43.170264850848675)),
 (0.518, (1369226.2668632688, 43.164744871516575)),
 (0.519, (1371225.1358790824, 43.15623425368964)),
 (0.52, (1371225.1358790824, 43.15623425368964)),
 (0.521, (1371225.1358790824, 43.15623425368964)),
 (0.522, (1371225.1358790824, 43.15623425368964)),
 (0.523, (1379684.1710197094, 43.12057509147576)),
 (0.524, (1381474.8082966425, 43.11310198538347)),
 (0.525, (1382983.9538093756, 43.10684997571817)),
 (0.526, (1382983.9538093756, 43.10684997571817)),
 (0.527, (1382983.9538093756, 43.10684997571817)),
 (0.528, (1382983.9538093756, 43.10684997571817)),
 (0.529, (1390178.0820348314, 43.077420885256956)),
 (0.53, (1392637.6176728825, 43.06748290358141)),
 (0.531, (1393580.7680392403, 43.063700851245)),
 (0.532, (1395724.8809742378, 43.05523250733681)),
 (0.533, (1395724.8809742378, 43.05523250733681)),
 (0.534, (1399209.3209597277, 43.04166621512637)),
 (0.535, (1401659.1241802094, 43.032224895914005)),
 (0.536, (1401659.1241802094, 43.032224895914005)),
 (0.537, (1401659.1241802094, 43.032224895914005)),
 (0.538, (1407576.2034714874, 43.00974246840563)),
 (0.539, (1408570.9840311836, 43.00602573388254)),
 (0.54, (1410183.332897933, 43.00005179299807)),
 (0.541, (1413018.9791654502, 42.98958751013153)),
 (0.542, (1414350.474571311, 42.98471479267136)),
 (0.543, (1416539.9722347488, 42.976742655459375)),
 (0.544, (1418769.6295302422, 42.96873708721449)),
 (0.545, (1418769.6295302422, 42.96873708721449)),
 (0.546, (1422103.331780041, 42.956970753208424)),
 (0.547, (1424095.9210660155, 42.94998358416385)),
 (0.548, (1425641.6978273685, 42.94462349600489)),
 (0.549, (1425641.6978273685, 42.94462349600489)),
 (0.55, (1429385.5304031502, 42.931763840189795)),
 (0.551, (1429385.5304031502, 42.931763840189795)),
 (0.552, (1429385.5304031502, 42.931763840189795)),
 (0.553, (1429385.5304031502, 42.931763840189795)),
 (0.554, (1429385.5304031502, 42.931763840189795)),
 (0.555, (1438727.3242142356, 42.900073390182754)),
 (0.556, (1440553.9174341452, 42.893959657063085)),
 (0.557, (1440553.9174341452, 42.893959657063085)),
 (0.558, (1443626.7962431898, 42.8837379280005)),
 (0.559, (1443626.7962431898, 42.8837379280005)),
 (0.56, (1443626.7962431898, 42.8837379280005)),
 (0.561, (1449682.5906128837, 42.863927347089884)),
 (0.562, (1449682.5906128837, 42.863927347089884)),
 (0.5630000000000001, (1453311.2794306476, 42.85217136287191)),
 (0.5640000000000001, (1455330.644484309, 42.84568444836188)),
 (0.5650000000000001, (1455330.644484309, 42.84568444836188)),
 (0.5660000000000001, (1457597.784555458, 42.83842963168615)),
 (0.5670000000000001, (1457597.784555458, 42.83842963168615)),
 (0.5680000000000001, (1457597.784555458, 42.83842963168615)),
 (0.5690000000000001, (1457597.784555458, 42.83842963168615)),
 (0.5700000000000001, (1457597.784555458, 42.83842963168615)),
 (0.5710000000000001, (1467881.957072041, 42.80593128438942)),
 (0.5720000000000001, (1470047.5313234595, 42.79914877141473)),
 (0.5730000000000001, (1471355.1205017003, 42.79507760741269)),
 (0.5740000000000001, (1473560.5992964958, 42.788234306404)),
 (0.5750000000000001, (1474193.7798087061, 42.78627961043839)),
 (0.5760000000000001, (1477525.6732188696, 42.776063627209965)),
 (0.577, (1478571.6090527945, 42.77287199343539)),
 (0.578, (1478571.6090527945, 42.77287199343539)),
 (0.579, (1482588.1332066883, 42.760672404295576)),
 (0.58, (1484449.3384963144, 42.75508867557045)),
 (0.581, (1484449.3384963144, 42.75508867557045)),
 (0.582, (1484449.3384963144, 42.75508867557045)),
 (0.583, (1489495.3457282702, 42.74014324052423)),
 (0.584, (1492086.0520159716, 42.73259212611171)),
 (0.585, (1492086.0520159716, 42.73259212611171)),
 (0.586, (1494417.5546012283, 42.72586929656484)),
 (0.587, (1496489.8826967315, 42.7199550120354)),
 (0.588, (1499091.326135834, 42.712598615773274)),
 (0.589, (1499091.326135834, 42.712598615773274)),
 (0.59, (1502356.9655342377, 42.70342024945375)),
 (0.591, (1504888.1505252013, 42.696377683037234)),
 (0.592, (1504888.1505252013, 42.696377683037234)),
 (0.593, (1504888.1505252013, 42.696377683037234)),
 (0.594, (1510226.1206418383, 42.68190693511459)),
 (0.595, (1511644.406202314, 42.67809243631294)),
 (0.596, (1511644.406202314, 42.67809243631294)),
 (0.597, (1514154.245903524, 42.67141668482116)),
 (0.598, (1516312.1208041508, 42.665750354003485)),
 (0.599, (1516312.1208041508, 42.665750354003485)),
 (0.6, (1516312.1208041508, 42.665750354003485)),
 (0.601, (1521283.1644316472, 42.65278934369559)),
 (0.602, (1521283.1644316472, 42.65278934369559)),
 (0.603, (1521283.1644316472, 42.65278934369559)),
 (0.604, (1528000.55015498, 42.63548013845723)),
 (0.605, (1529954.6982695658, 42.63050807786182)),
 (0.606, (1531489.568651363, 42.626701930636486)),
 (0.607, (1532607.031279209, 42.623948443914564)),
 (0.608, (1535031.7208215182, 42.61801482605009)),
 (0.609, (1535031.7208215182, 42.61801482605009)),
 (0.61, (1537767.604763809, 42.61141040249598)),
 (0.611, (1540508.40605797, 42.604878302856456)),
 (0.612, (1540508.40605797, 42.604878302856456)),
 (0.613, (1544465.101935842, 42.59551300004948)),
 (0.614, (1546129.3594892966, 42.59163146959967)),
 (0.615, (1547217.6763131209, 42.589121661685795)),
 (0.616, (1547217.6763131209, 42.589121661685795)),
 (0.617, (1547217.6763131209, 42.589121661685795)),
 (0.618, (1547217.6763131209, 42.589121661685795)),
 (0.619, (1547217.6763131209, 42.589121661685795)),
 (0.62, (1556131.8537093808, 42.56883223033582)),
 (0.621, (1556131.8537093808, 42.56883223033582)),
 (0.622, (1556131.8537093808, 42.56883223033582)),
 (0.623, (1556131.8537093808, 42.56883223033582)),
 (0.624, (1556131.8537093808, 42.56883223033582)),
 (0.625, (1556131.8537093808, 42.56883223033582)),
 (0.626, (1556131.8537093808, 42.56883223033582)),
 (0.627, (1569155.3380843708, 42.53987870805724)),
 (0.628, (1570973.7676858357, 42.53591159966429)),
 (0.629, (1572268.025468364, 42.5330995075112)),
 (0.63, (1572268.025468364, 42.5330995075112)),
 (0.631, (1572268.025468364, 42.5330995075112)),
 (0.632, (1577819.5551494425, 42.52122845245574)),
 (0.633, (1579207.9675316876, 42.51829998991907)),
 (0.634, (1579207.9675316876, 42.51829998991907)),
 (0.635, (1579207.9675316876, 42.51829998991907)),
 (0.636, (1579207.9675316876, 42.51829998991907)),
 (0.637, (1586054.6946943023, 42.50398151178249)),
 (0.638, (1588310.9521224909, 42.499297114068604)),
 (0.639, (1590811.513966618, 42.4941360326921)),
 (0.64, (1592355.9299020236, 42.49099505858163)),
 (0.641, (1593520.3639331204, 42.48864838581197)),
 (0.642, (1596283.7490394549, 42.48314935315228)),
 (0.643, (1596283.7490394549, 42.48314935315228)),
 (0.644, (1599118.192173848, 42.477589287323795)),
 (0.645, (1599118.192173848, 42.477589287323795)),
 (0.646, (1599118.192173848, 42.477589287323795)),
 (0.647, (1605401.5033365383, 42.46537704140558)),
 (0.648, (1605401.5033365383, 42.46537704140558)),
 (0.649, (1608095.9163783656, 42.46020305128399)),
 (0.65, (1608095.9163783656, 42.46020305128399)),
 (0.651, (1608095.9163783656, 42.46020305128399)),
 (0.652, (1613386.584077212, 42.4502984543921)),
 (0.653, (1613386.584077212, 42.4502984543921)),
 (0.654, (1613386.584077212, 42.4502984543921)),
 (0.655, (1613386.584077212, 42.4502984543921)),
 (0.656, (1613386.584077212, 42.4502984543921)),
 (0.657, (1621424.765186395, 42.43539594106051)),
 (0.658, (1623643.8370462516, 42.43131582900795)),
 (0.659, (1623643.8370462516, 42.43131582900795)),
 (0.66, (1627527.3887331164, 42.42427426801734)),
 (0.661, (1629327.8551135068, 42.42104620273736)),
 (0.662, (1631776.9930122467, 42.4167069548855)),
 (0.663, (1631776.9930122467, 42.4167069548855)),
 (0.664, (1631776.9930122467, 42.4167069548855)),
 (0.665, (1636275.3971363776, 42.40879979372882)),
 (0.666, (1639006.0228085304, 42.40413740995934)),
 (0.667, (1640703.0528782867, 42.40125933598711)),
 (0.668, (1642513.4217077591, 42.3982063220022)),
 (0.669, (1642513.4217077591, 42.3982063220022)),
 (0.67, (1642513.4217077591, 42.3982063220022)),
 (0.671, (1647247.6287427414, 42.39027797395335)),
 (0.672, (1649074.0824638307, 42.38727193226452)),
 (0.673, (1651435.4266256928, 42.38345125399764)),
 (0.674, (1646936.0852463106, 42.46163571202437)),
 (0.675, (1596803.4519518015, 42.45672631658782)),
 (0.676, (1599790.8715948632, 42.451761598726776)),
 (0.677, (1600781.3036055283, 42.450135931485825)),
 (0.678, (1600781.3036055283, 42.450135931485825)),
 (0.679, (1600781.3036055283, 42.450135931485825)),
 (0.68, (1600781.3036055283, 42.450135931485825)),
 (0.681, (1600781.3036055283, 42.450135931485825)),
 (0.682, (1600781.3036055283, 42.450135931485825)),
 (0.683, (1600781.3036055283, 42.450135931485825)),
 (0.684, (1600781.3036055283, 42.450135931485825)),
 (0.685, (1600781.3036055283, 42.450135931485825)),
 (0.686, (1600781.3036055283, 42.450135931485825)),
 (0.687, (1600781.3036055283, 42.450135931485825)),
 (0.6880000000000001, (1600781.3036055283, 42.450135931485825)),
 (0.6890000000000001, (1600781.3036055283, 42.450135931485825)),
 (0.6900000000000001, (1600781.3036055283, 42.450135931485825)),
 (0.6910000000000001, (1600781.3036055283, 42.450135931485825)),
 (0.6920000000000001, (1600781.3036055283, 42.450135931485825)),
 (0.6930000000000001, (1600781.3036055283, 42.450135931485825)),
 (0.6940000000000001, (1600781.3036055283, 42.450135931485825)),
 (0.6950000000000001, (1600781.3036055283, 42.450135931485825)),
 (0.6960000000000001, (1600781.3036055283, 42.450135931485825)),
 (0.6970000000000001, (1600781.3036055283, 42.450135931485825)),
 (0.6980000000000001, (1600781.3036055283, 42.450135931485825)),
 (0.6990000000000001, (1600781.3036055283, 42.450135931485825)),
 (0.7000000000000001, (1600781.3036055283, 42.450135931485825)),
 (0.7010000000000001, (1600781.3036055283, 42.450135931485825)),
 (0.7020000000000001, (1600781.3036055283, 42.450135931485825)),
 (0.7030000000000001, (1600781.3036055283, 42.450135931485825)),
 (0.704, (1600781.3036055283, 42.450135931485825)),
 (0.705, (1600781.3036055283, 42.450135931485825)),
 (0.706, (1597514.3036055283, 42.20825302448651)),
 (0.707, (1597514.3036055283, 42.20825302448651)),
 (0.708, (1597514.3036055283, 42.20825302448651)),
 (0.709, (1597514.3036055283, 42.20825302448651)),
 (0.71, (1597514.3036055283, 42.20825302448651)),
 (0.711, (1597514.3036055283, 42.20825302448651)),
 (0.712, (1597514.3036055283, 42.20825302448651)),
 (0.713, (1661346.4773899356, 42.1171811454383)),
 (0.714, (1658936.367716731, 42.19183232321753)),
 (0.715, (1645166.127716731, 42.29412572907674)),
 (0.716, (1648303.8016055124, 42.28954262265229)),
 (0.717, (1650044.3829549288, 42.28708710394712)),
 (0.718, (1650928.2129560313, 42.285869266164624)),
 (0.719, (1654002.3016276627, 42.28166071596513)),
 (0.72, (1644744.6249621569, 40.71643451371214)),
 (0.721, (1646642.6656229037, 40.715705551533524)),
 (0.722, (1647353.044084335, 40.71543864277567)),
 (0.723, (1650152.881095927, 40.71440298321249)),
 (0.724, (1650621.739738017, 40.714242397562195)),
 (0.725, (1653015.0989600557, 40.71345367576092)),
 (0.726, (1653015.0989600557, 40.71345367576092)),
 (0.727, (1655901.623090563, 40.71253776382341)),
 (0.728, (1655901.623090563, 40.71253776382341)),
 (0.729, (1660441.162711699, 40.71119680175664)),
 (0.73, (1662132.6269994548, 40.7107660569616)),
 (0.731, (1662132.6269994548, 40.7107660569616)),
 (0.732, (1665589.7387862017, 40.709913446061165)),
 (0.733, (1667339.4993075284, 40.70949949576115)),
 (0.734, (1669169.1785453279, 40.70909818635682)),
 (0.735, (1670529.810126635, 40.70880610735869)),
 (0.736, (1670529.810126635, 40.70880610735869)),
 (0.737, (1670529.810126635, 40.70880610735869)),
 (0.738, (1675310.0249780242, 40.707838279338326)),
 (0.739, (1676688.63114592, 40.70756642913803)),
 (0.74, (1678525.395965629, 40.7072159299589)),
 (0.741, (1678525.395965629, 40.7072159299589)),
 (0.742, (1682747.6938063533, 40.70643498813516)),
 (0.743, (1684254.285889074, 40.70617898799609)),
 (0.744, (1686554.0413080982, 40.70580458695191)),
 (0.745, (1686554.0413080982, 40.70580458695191)),
 (0.746, (1689521.9884616295, 40.7054266326454)),
 (0.747, (1691585.3810323426, 40.705185956203955)),
 (0.748, (1691585.3810323426, 40.705185956203955)),
 (0.749, (1694215.876922184, 40.70489601899619)),
 (0.75, (1694215.876922184, 40.70489601899619)),
 (0.751, (1698686.3964483014, 40.704448633109095)),
 (0.752, (1699402.7382996231, 40.70437717842187)),
 (0.753, (1699402.7382996231, 40.70437717842187)),
 (0.754, (1702828.5672207854, 40.70406598825349)),
 (0.755, (1704281.0093980301, 40.70393793106664)),
 (0.756, (1706256.0487849286, 40.70378028452452)),
 (0.757, (1708959.3634961592, 40.70365248965767)),
 (0.758, (1708959.3634961592, 40.70365248965767)),
 (0.759, (1708959.3634961592, 40.70365248965767)),
 (0.76, (1713394.047050823, 40.70347049496915)),
 (0.761, (1715307.131979079, 40.70345814097489)),
 (0.762, (1716594.8555885812, 40.70347088202108)),
 (0.763, (1718227.3410025674, 40.70349885131586)),
 (0.764, (1716133.3410025674, 40.55907762634134)),
 (0.765, (1719537.6446898018, 40.55944496168652)),
 (0.766, (1722365.5342381606, 40.559840068064574)),
 (0.767, (1722365.5342381606, 40.559840068064574)),
 (0.768, (1716897.3742381604, 40.505705368205604)),
 (0.769, (1716897.3742381604, 40.505705368205604)),
 (0.77, (1723572.5322490223, 40.50689139352863)),
 (0.771, (1725194.647061754, 40.50723824293909)),
 (0.772, (1727239.0620245372, 40.50779449400937)),
 (0.773, (1725108.1020245373, 39.67365096254044)),
 (0.774, (1725108.1020245373, 39.67365096254044)),
 (0.775, (1725108.1020245373, 39.67365096254044)),
 (0.776, (1731795.0810302733, 39.67879641073012)),
 (0.777, (1731795.0810302733, 39.67879641073012)),
 (0.778, (1731795.0810302733, 39.67879641073012)),
 (0.779, (1731795.0810302733, 39.67879641073012)),
 (0.78, (1738662.9980616034, 39.684200292089315)),
 (0.781, (1738662.9980616034, 39.684200292089315)),
 (0.782, (1738662.9980616034, 39.684200292089315)),
 (0.783, (1738662.9980616034, 39.684200292089315)),
 (0.784, (1745300.3186967513, 39.68952017967793)),
 (0.785, (1747109.219937373, 39.69099813468368)),
 (0.786, (1748405.3441451762, 39.69208365498057)),
 (0.787, (1750580.270896708, 39.693922843051055)),
 (0.788, (1752296.3701315718, 39.69542084430647)),
 (0.789, (1753661.083809007, 39.69661420480114)),
 (0.79, (1754427.6388631272, 39.69728516127052)),
 (0.791, (1757386.4464879336, 39.69992018752057)),
 (0.792, (1758530.7287586448, 39.70094259381216)),
 (0.793, (1758530.7287586448, 39.70094259381216)),
 (0.794, (1761154.3537125774, 39.70331347222855)),
 (0.795, (1760126.1937125775, 39.623622361900786)),
 (0.796, (1764689.9425301438, 39.62797705876281)),
 (0.797, (1764346.7716519544, 39.25399652144479)),
 (0.798, (1764346.7716519544, 39.25399652144479)),
 (0.799, (1768071.061551408, 39.258466710385235)),
 (0.8, (1769985.9629860106, 39.260780160914685)),
 (0.801, (1770425.5686484114, 39.261316270282684)),
 (0.802, (1773283.933414241, 39.2648652000591)),
 (0.803, (1775261.5170704531, 39.26736515233705)),
 (0.804, (1775261.5170704531, 39.26736515233705)),
 (0.805, (1775261.5170704531, 39.26736515233705)),
 (0.806, (1779991.0873030403, 39.2733453020131)),
 (0.807, (1779991.0873030403, 39.2733453020131)),
 (0.808, (1775803.3873030401, 39.310605732826104)),
 (0.809, (1763024.6624059796, 36.6453410683926)),
 (0.81, (1763024.6624059796, 36.6453410683926)),
 (0.811, (1765714.1884242543, 36.65287082355967)),
 (0.812, (1768180.751185677, 36.659788835880526)),
 (0.8130000000000001, (1769733.8363652243, 36.66416090940172)),
 (0.8140000000000001, (1770129.516755981, 36.6652754137424)),
 (0.8150000000000001, (1770129.516755981, 36.6652754137424)),
 (0.8160000000000001, (1770129.516755981, 36.6652754137424)),
 (0.8170000000000001, (1770129.516755981, 36.6652754137424)),
 (0.8180000000000001, (1770129.516755981, 36.6652754137424)),
 (0.8190000000000001, (1778861.9678809023, 36.68991435721737)),
 (0.8200000000000001, (1778861.9678809023, 36.68991435721737)),
 (0.8210000000000001, (1778861.9678809023, 36.68991435721737)),
 (0.8220000000000001, (1778861.9678809023, 36.68991435721737)),
 (0.8230000000000001, (1778861.9678809023, 36.68991435721737)),
 (0.8240000000000001, (1787099.4219535603, 36.71296000364977)),
 (0.8250000000000001, (1789950.2297993782, 36.72094460144245)),
 (0.8260000000000001, (1791334.1586912917, 36.724817165135235)),
 (0.8270000000000001, (1792705.260634783, 36.728679885504036)),
 (0.8280000000000001, (1794687.9799656023, 36.734269771511165)),
 (0.8290000000000001, (1796401.6225554068, 36.7391056365538)),
 (0.8300000000000001, (1796401.6225554068, 36.7391056365538)),
 (0.8310000000000001, (1799445.5261900988, 36.74767489458937)),
 (0.8320000000000001, (1800375.3827280044, 36.75028784883137)),
 (0.833, (1800375.3827280044, 36.75028784883137)),
 (0.834, (1800375.3827280044, 36.75028784883137)),
 (0.835, (1800375.3827280044, 36.75028784883137)),
 (0.836, (1800375.3827280044, 36.75028784883137)),
 (0.837, (1800375.3827280044, 36.75028784883137)),
 (0.838, (1810647.0467280042, 36.77900399388079)),
 (0.839, (1812703.6892941808, 36.784724374726736)),
 (0.84, (1814976.1802845176, 36.791057503851505)),
 (0.841, (1814976.1802845176, 36.791057503851505)),
 (0.842, (1814976.1802845176, 36.791057503851505)),
 (0.843, (1814976.1802845176, 36.791057503851505)),
 (0.844, (1820457.991637114, 36.8064294368114)),
 (0.845, (1820457.991637114, 36.8064294368114)),
 (0.846, (1824890.7682174442, 36.81882838226799)),
 (0.847, (1824890.7682174442, 36.81882838226799)),
 (0.848, (1828532.6275552162, 36.82900528819303)),
 (0.849, (1830014.0049747524, 36.83314743439176)),
 (0.85, (1831515.7481069143, 36.83734330212511)),
 (0.851, (1832234.3585115648, 36.83935053679468)),
 (0.852, (1832234.3585115648, 36.83935053679468)),
 (0.853, (1832234.3585115648, 36.83935053679468)),
 (0.854, (1837723.499598943, 36.85463449820539)),
 (0.855, (1837723.499598943, 36.85463449820539)),
 (0.856, (1837723.499598943, 36.85463449820539)),
 (0.857, (1842662.2640486788, 36.86839934015767)),
 (0.858, (1842662.2640486788, 36.86839934015767)),
 (0.859, (1842662.2640486788, 36.86839934015767)),
 (0.86, (1842662.2640486788, 36.86839934015767)),
 (0.861, (1842662.2640486788, 36.86839934015767)),
 (0.862, (1851942.8262612366, 36.894125223002156)),
 (0.863, (1852135.6767601073, 36.894657900393334)),
 (0.864, (1852135.6767601073, 36.894657900393334)),
 (0.865, (1855595.8370782284, 36.904216439234766)),
 (0.866, (1855595.8370782284, 36.904216439234766)),
 (0.867, (1860237.0736111205, 36.91714071901851)),
 (0.868, (1860237.0736111205, 36.91714071901851)),
 (0.869, (1860237.0736111205, 36.91714071901851)),
 (0.87, (1864108.3232707994, 36.92787904377297)),
 (0.871, (1864108.3232707994, 36.92787904377297)),
 (0.872, (1868456.9050696285, 36.93989816318943)),
 (0.873, (1868456.9050696285, 36.93989816318943)),
 (0.874, (1868456.9050696285, 36.93989816318943)),
 (0.875, (1873296.6461996494, 36.953225025194975)),
 (0.876, (1875233.8129750078, 36.958543090014985)),
 (0.877, (1875981.297129798, 36.96059545377948)),
 (0.878, (1875981.297129798, 36.96059545377948)),
 (0.879, (1880047.6159740987, 36.971774175475545)),
 (0.88, (1881562.9510791958, 36.97600346467724)),
 (0.881, (1882826.6764298012, 36.97953703965471)),
 (0.882, (1884706.2278034957, 36.98478987502741)),
 (0.883, (1885736.1269511934, 36.98766655306129)),
 (0.884, (1887776.5309306253, 36.99335701489712)),
 (0.885, (1887776.5309306253, 36.99335701489712)),
 (0.886, (1888923.757798214, 34.66762236519408)),
 (0.887, (1890109.71013542, 34.672420872793325)),
 (0.888, (1890109.71013542, 34.672420872793325)),
 (0.889, (1890109.71013542, 34.672420872793325)),
 (0.89, (1890109.71013542, 34.672420872793325)),
 (0.891, (1890109.71013542, 34.672420872793325)),
 (0.892, (1890109.71013542, 34.672420872793325)),
 (0.893, (1890109.71013542, 34.672420872793325)),
 (0.894, (1890109.71013542, 34.672420872793325)),
 (0.895, (1903501.1583506772, 34.72634320775449)),
 (0.896, (1903501.1583506772, 34.72634320775449)),
 (0.897, (1906667.120098088, 34.73901515274373)),
 (0.898, (1906667.120098088, 34.73901515274373)),
 (0.899, (1909499.4747170806, 34.7503570799983)),
 (0.9, (1911710.3133978467, 34.75920835602772)),
 (0.901, (1911710.3133978467, 34.75920835602772)),
 (0.902, (1911710.3133978467, 34.75920835602772)),
 (0.903, (1911710.3133978467, 34.75920835602772)),
 (0.904, (1918685.6745838332, 34.787058327585925)),
 (0.905, (1919892.7406931673, 34.79190090697813)),
 (0.906, (1921779.8193519516, 34.799503114307655)),
 (0.907, (1922058.9918212378, 34.800630868268556)),
 (0.908, (1922058.9918212378, 34.800630868268556)),
 (0.909, (1926215.3556158002, 34.81744854195251)),
 (0.91, (1928403.22378676, 34.826276161846344)),
 (0.911, (1929527.4684798832, 34.83081403729165)),
 (0.912, (1931133.8798334722, 34.83729547615242)),
 (0.913, (1932816.6012057953, 34.844080641229986)),
 (0.914, (1932816.6012057953, 34.844080641229986)),
 (0.915, (1932816.6012057953, 34.844080641229986)),
 (0.916, (1938453.6191123489, 34.86673528351463)),
 (0.917, (1939338.9517181413, 34.8702929740858)),
 (0.918, (1941606.9948264377, 34.87940643946802)),
 (0.919, (1941803.0662233133, 34.88019517015477)),
 (0.92, (1944889.9160013136, 34.89261192770758)),
 (0.921, (1945871.9564044562, 34.89655438732529)),
 (0.922, (1945871.9564044562, 34.89655438732529)),
 (0.923, (1948372.1057341264, 34.90657756583472)),
 (0.924, (1951479.4615389556, 34.91900330190154)),
 (0.925, (1951479.4615389556, 34.91900330190154)),
 (0.926, (1951479.4615389556, 34.91900330190154)),
 (0.927, (1955338.5714012883, 34.93443831587016)),
 (0.928, (1957912.6680791155, 34.94471323919435)),
 (0.929, (1959411.7011460653, 34.95070896558495)),
 (0.93, (1959411.7011460653, 34.95070896558495)),
 (0.931, (1959411.7011460653, 34.95070896558495)),
 (0.932, (1963852.4948032387, 34.96843387235004)),
 (0.933, (1965400.5284358142, 34.97461571386869)),
 (0.934, (1967023.3363093412, 34.98110335771664)),
 (0.935, (1969058.006592242, 34.989324934657056)),
 (0.936, (1970945.3545428757, 34.996961168495766)),
 (0.937, (1972340.9155269475, 35.002629905345074)),
 (0.9380000000000001, (1973403.669157547, 35.00696090561922)),
 (0.9390000000000001, (1973403.669157547, 35.00696090561922)),
 (0.9400000000000001, (1976253.3243626468, 35.01860190130032)),
 (0.9410000000000001, (1979076.017567865, 35.03010191082997)),
 (0.9420000000000001, (1979832.417567865, 35.03317854622332)),
 (0.9430000000000001, (1979832.417567865, 35.03317854622332)),
 (0.9440000000000001, (1983179.219007039, 35.046783881052946)),
 (0.9450000000000001, (1983179.219007039, 35.046783881052946)),
 (0.9460000000000001, (1986625.9491780247, 35.06076035794786)),
 (0.9470000000000001, (1988100.6967706233, 35.066728732834584)),
 (0.9480000000000001, (1988100.6967706233, 35.066728732834584)),
 (0.9490000000000001, (1992210.0356570226, 35.08335649798272)),
 (0.9500000000000001, (1992656.1225223537, 35.085166517136756)),
 (0.9510000000000001, (1993982.6245453886, 35.08984287267893)),
 (0.9520000000000001, (1991219.6438592963, 35.070661605078186)),
 (0.9530000000000001, (1991219.6438592963, 35.070661605078186)),
 (0.9540000000000001, (1990120.2038592964, 35.00450208191391)),
 (0.9550000000000001, (1995569.2289316377, 35.026942873248075)),
 (0.9560000000000001, (1993006.6534266125, 34.07481011099974)),
 (0.9570000000000001, (1994077.7219043418, 34.07973163457018)),
 (0.9580000000000001, (1996344.4802670411, 34.09013744578911)),
 (0.9590000000000001, (1996344.4802670411, 34.09013744578911)),
 (0.96, (1998247.4800384867, 34.09885676821181)),
 (0.961, (2000844.9055861712, 34.11074420674395)),
 (0.962, (2001760.3819523714, 34.11493664528704)),
 (0.963, (2004013.9959245203, 34.125259836110615)),
 (0.964, (2004625.4759434313, 34.128060315732476)),
 (0.965, (2007378.8800907831, 34.14064947501463)),
 (0.966, (2008654.8347714012, 34.14647274480419)),
 (0.967, (2010770.7938715508, 34.1561266898063)),
 (0.968, (2012577.0010369993, 34.164383459511356)),
 (0.969, (2013847.8634748417, 34.17019617159261)),
 (0.97, (2015630.7282606578, 34.17835952480938)),
 (0.971, (2017313.0928216837, 34.186060916855794)),
 (0.972, (2018813.6293186292, 34.19292657343878)),
 (0.973, (2020706.499873178, 34.20160340128224)),
 (0.974, (2020706.499873178, 34.20160340128224)),
 (0.975, (2023437.2064653006, 34.21410310235987)),
 (0.976, (2024848.8874757807, 34.22055267034169)),
 (0.977, (2026407.0671254068, 34.22766623810515)),
 (0.978, (2028756.7199919492, 34.23839413086611)),
 (0.979, (2028756.7199919492, 34.23839413086611)),
 (0.98, (2028756.7199919492, 34.23839413086611)),
 (0.981, (2033288.7490129452, 34.25905252523253)),
 (0.982, (2034224.0850522218, 34.26330945258805)),
 (0.983, (2036264.8484946967, 34.27260845452638)),
 (0.984, (2036853.5869306035, 34.27529108400344)),
 (0.985, (2039837.4389131628, 34.28887613024939)),
 (0.986, (2033571.818913163, 34.16509082720591)),
 (0.987, (2035863.1859161742, 34.1756699952195)),
 (0.988, (2037827.1135792038, 34.18473400046884)),
 (0.989, (2039600.6913852254, 34.19290503419731)),
 (0.99, (2040269.0443053609, 34.195992266346465)),
 (0.991, (2043380.1520602964, 34.21035401725523)),
 (0.992, (2043380.1520602964, 34.21035401725523)),
 (0.993, (2042852.781002447, 33.997195544076384)),
 (0.994, (2045439.5961530032, 34.00938838201095)),
 (0.995, (2045840.7266019671, 34.01128039611291)),
 (0.996, (2045840.7266019671, 34.01128039611291)),
 (0.997, (2050142.421807609, 34.031545722549545)),
 (0.998, (2050142.421807609, 34.031545722549545)),
 (0.999, (2053765.617620457, 34.04859254038655)),
 ...]


# Extracting total CO2 avoided and average cost per tCO2 for each allocation
total_co2_avoided2 = [result[0] for _, result in results_with_allocations2]
avg_cost_per_tco22 = [result[1] for _, result in results_with_allocations2]

# Plotting the results
fig, ax12 = plt.subplots()

# Bar graph for total CO2 avoided
ax12.bar(allocations, total_co2_avoided2, width=0.03, color='b', align='center')
ax12.set_xlabel('Allocation to Scenario 1')
ax12.set_ylabel('Total CO2 Avoided (tCO2)', color='b')
ax12.tick_params('y', colors='b')

# Creating a second y-axis for average cost per tCO2
ax13 = ax12.twinx()
ax13.plot(allocations, avg_cost_per_tco22, 'r-')
ax13.set_ylabel('Average Cost per tCO2 (USD)', color='r')
ax13.tick_params('y', colors='r')

plt.title('Total CO2 Avoided and Average Cost per tCO2 by Allocation')
plt.show()


dm['tax_offer'] = dm['acq_cost_pred_usd'] * 0.5
dm['tax_offer'] = dm['tax_offer'].clip(upper=75000)
dm['flatrate'] = dm['p_f_2010']/100 * dm['ha'] * 1210
i_scenario_for = dm['tax_offer'] < dm['flatrate']
selected_parcels_scenario_for = dm[i_scenario_for]
dm['carbonrate']= dm['tco2_loss_pred_100yr'] *13.1
i_scenario_carbonrate = dm['tax_offer'] < dm['carbonrate']


def simulate_allocation(budget_allocation, dm, BUDGET):
    """
    Simulates the allocation of budget between the two scenarios and calculates the total CO2 avoided 
    and average cost per tCO2 for the combined allocation.
    """
    # Allocating budget according to the given ratio
    budget_scenario_for = BUDGET * budget_allocation
    budget_scenario_car = BUDGET * (1 - budget_allocation)

    # Scenario 1
    i_scenario_for = dm['flatrate'].sum().lt(budget_scenario_for)
    selected_scenario_for = dm[i_scenario_for]

    # Scenario 2
    i_scenario_car = i_scenario_carbonrate['carbonrate'].sum().lt(budget_scenario_car)
    selected_scenario_car = dm[i_scenario_car]

    # Combined results
    combined_selected = pd.concat([selected_scenario_for, selected_scenario_car])
    total_co2_avoided = combined_selected['tco2_loss_pred_100yr'].sum()
    total_budget_used = combined_selected['ce_cost_pred_usd'].sum()
    avg_cost_per_tco2 = total_budget_used / total_co2_avoided if total_co2_avoided > 0 else np.inf

    return total_co2_avoided, avg_cost_per_tco2

allocations = np.arange(0, 1.05, 0.05)  # Allocating from 0% to 100% in 5% increments for Scenario 1
results3 = [simulate_allocation(allocation, dm, BUDGET) for allocation in allocations]

# Pairing results with allocations
results_with_allocations3 = list(zip(allocations, results))
results_with_allocations3

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
File ~/anaconda3/lib/python3.11/site-packages/pandas/core/indexes/base.py:3653, in Index.get_loc(self, key)
   3652 try:
-> 3653     return self._engine.get_loc(casted_key)
   3654 except KeyError as err:

File ~/anaconda3/lib/python3.11/site-packages/pandas/_libs/index.pyx:147, in pandas._libs.index.IndexEngine.get_loc()

File ~/anaconda3/lib/python3.11/site-packages/pandas/_libs/index.pyx:176, in pandas._libs.index.IndexEngine.get_loc()

File pandas/_libs/hashtable_class_helper.pxi:7080, in pandas._libs.hashtable.PyObjectHashTable.get_item()

File pandas/_libs/hashtable_class_helper.pxi:7088, in pandas._libs.hashtable.PyObjectHashTable.get_item()

KeyError: 'flatrate'

The above exception was the direct cause of the following exception:

KeyError                                  Traceback (most recent call last)
Cell In[159], line 27
     24     return total_co2_avoided, avg_cost_per_tco2
     26 allocations = np.arange(0, 1.05, 0.05)  # Allocating from 0% to 100% in 5% increments for Scenario 1
---> 27 results3 = [simulate_allocation(allocation, dm, BUDGET) for allocation in allocations]
     29 # Pairing results with allocations
     30 results_with_allocations3 = list(zip(allocations, results))

Cell In[159], line 27, in <listcomp>(.0)
     24     return total_co2_avoided, avg_cost_per_tco2
     26 allocations = np.arange(0, 1.05, 0.05)  # Allocating from 0% to 100% in 5% increments for Scenario 1
---> 27 results3 = [simulate_allocation(allocation, dm, BUDGET) for allocation in allocations]
     29 # Pairing results with allocations
     30 results_with_allocations3 = list(zip(allocations, results))

Cell In[159], line 11, in simulate_allocation(budget_allocation, dm, BUDGET)
      8 budget_scenario_car = BUDGET * (1 - budget_allocation)
     10 # Scenario 1
---> 11 i_scenario_for_sum = i_scenario_for['flatrate'].sum().lt(budget_scenario_for)
     12 selected_scenario_for = dm[i_scenario_for_sum]
     14 # Scenario 2

File ~/anaconda3/lib/python3.11/site-packages/pandas/core/series.py:1007, in Series.__getitem__(self, key)
   1004     return self._values[key]
   1006 elif key_is_scalar:
-> 1007     return self._get_value(key)
   1009 if is_hashable(key):
   1010     # Otherwise index.get_value will raise InvalidIndexError
   1011     try:
   1012         # For labels that don't resolve as scalars like tuples and frozensets

File ~/anaconda3/lib/python3.11/site-packages/pandas/core/series.py:1116, in Series._get_value(self, label, takeable)
   1113     return self._values[label]
   1115 # Similar to Index.get_value, but we do not fall back to positional
-> 1116 loc = self.index.get_loc(label)
   1118 if is_integer(loc):
   1119     return self._values[loc]

File ~/anaconda3/lib/python3.11/site-packages/pandas/core/indexes/base.py:3655, in Index.get_loc(self, key)
   3653     return self._engine.get_loc(casted_key)
   3654 except KeyError as err:
-> 3655     raise KeyError(key) from err
   3656 except TypeError:
   3657     # If we have a listlike key, _check_indexing_error will raise
   3658     #  InvalidIndexError. Otherwise we fall through and re-raise
   3659     #  the TypeError.
   3660     self._check_indexing_error(key)

KeyError: 'flatrate'

Dep. Variable:	ls_price_per_ha_log	R-squared:	0.455
Model:	OLS	Adj. R-squared:	0.455
Method:	Least Squares	F-statistic:	1687.
Date:	Thu, 14 Dec 2023	Prob (F-statistic):	0.00
Time:	17:18:24	Log-Likelihood:	-78212.
No. Observations:	44424	AIC:	1.565e+05
Df Residuals:	44401	BIC:	1.567e+05
Df Model:	22
Covariance Type:	nonrobust

Omnibus:	7448.307	Durbin-Watson:	1.958
Prob(Omnibus):	0.000	Jarque-Bera (JB):	34364.698
Skew:	-0.752	Prob(JB):	0.00
Kurtosis:	7.038	Cond. No.	4.82e+07

Dep. Variable:	ls_price_per_ha_log	R-squared:	0.462
Model:	OLS	Adj. R-squared:	0.461
Method:	Least Squares	F-statistic:	1002.
Date:	Thu, 14 Dec 2023	Prob (F-statistic):	0.00
Time:	17:18:31	Log-Likelihood:	-77948.
No. Observations:	44424	AIC:	1.560e+05
Df Residuals:	44385	BIC:	1.563e+05
Df Model:	38
Covariance Type:	nonrobust

Omnibus:	7598.673	Durbin-Watson:	1.954
Prob(Omnibus):	0.000	Jarque-Bera (JB):	37366.196
Skew:	-0.749	Prob(JB):	0.00
Kurtosis:	7.236	Cond. No.	2.16e+06

Dep. Variable:	ls_price_per_ha_log	R-squared:	0.537
Model:	OLS	Adj. R-squared:	0.534
Method:	Least Squares	F-statistic:	190.2
Date:	Thu, 14 Dec 2023	Prob (F-statistic):	0.00
Time:	17:18:33	Log-Likelihood:	-74613.
No. Observations:	44424	AIC:	1.498e+05
Df Residuals:	44154	BIC:	1.521e+05
Df Model:	269
Covariance Type:	nonrobust

	coef	std err	t	P>\|t\|	[0.025	0.975]
const	-52.9732	2.823	-18.764	0.000	-58.507	-47.440
slope	-0.0266	0.002	-10.791	0.000	-0.031	-0.022
travel	-0.0006	0.000	-6.079	0.000	-0.001	-0.000
p_wet	-0.0050	0.000	-24.097	0.000	-0.005	-0.005
coast_2500	0.0298	0.001	29.846	0.000	0.028	0.032
river_frontage	-2.028e-05	7.75e-05	-0.262	0.794	-0.000	0.000
lake_frontage	-0.0007	5.09e-05	-12.880	0.000	-0.001	-0.001
pop_dens_tract_2012-2016	-3.011e-05	1.99e-05	-1.511	0.131	-6.92e-05	8.94e-06
hh_inc_avg_tract_2012-2016	1.322e-05	1.88e-07	70.248	0.000	1.29e-05	1.36e-05
p_f_2000	-0.0059	0.000	-16.005	0.000	-0.007	-0.005
p_prot_2000_500	-2.324e-05	0.001	-0.046	0.964	-0.001	0.001
p_d_2000_200	0.0124	0.001	17.975	0.000	0.011	0.014
p_d_2000_1000	-0.0002	0.001	-0.250	0.803	-0.002	0.002
p_d_2000_5000	0.0154	0.001	16.528	0.000	0.014	0.017
p_dev_open	0.0036	0.001	6.977	0.000	0.003	0.005
p_dev_low	0.0075	0.001	9.104	0.000	0.006	0.009
p_dev_medium	0.0150	0.001	20.060	0.000	0.013	0.016
p_dev_high	0.0045	0.001	5.636	0.000	0.003	0.006
p_crops	-0.0038	0.001	-4.302	0.000	-0.006	-0.002
p_pasture	-0.0002	0.001	-0.382	0.702	-0.001	0.001
p_grassland	0.0024	0.001	2.224	0.026	0.000	0.004
bld_area_per_ha	1.86e-05	6.64e-07	28.024	0.000	1.73e-05	1.99e-05
ls_year	0.0315	0.001	22.385	0.000	0.029	0.034

Omnibus:	8803.431	Durbin-Watson:	1.979
Prob(Omnibus):	0.000	Jarque-Bera (JB):	59018.463
Skew:	-0.788	Prob(JB):	0.00
Kurtosis:	8.422	Cond. No.	2.31e+15

	fips	csd_id	ha	slope	travel	p_wet	coast_2500	river_frontage	lake_frontage	lake_importance	...	hh_inc_med_tract_2012-2016	hh_inc_avg_bg_2012-2016	hh_inc_avg_tract_2012-2016	p_f_change	ls_price_per_ha_log	ls_year	p_e_ls	bld_area_per_ha	ls_pred	fc_pred
pid
F_803174_3027860	25009	16250	1.138	3.269	8.000	0.0	0.000	0.0	0.0	0.000	...	98264.0	107637.264151	111868.256722	0	NaN	NaN	0.0	3124.780316	12.870125	-2.555327
M_238118_942044	25009	27620	3.052	6.315	14.000	10.0	0.000	0.0	0.0	0.000	...	96512.0	127585.099685	115417.941953	0	9.598621	1991.0	0.0	1471.821756	12.174028	-4.203807
M_240445_943056	25009	27620	2.165	2.538	20.000	48.0	0.000	0.0	0.0	0.000	...	96512.0	127585.099685	115417.941953	0	NaN	NaN	0.0	3155.196305	12.654431	-1.659536
M_245474_951109	25009	77150	2.162	1.784	19.000	21.0	7.000	0.0	0.0	0.000	...	138947.0	168559.326425	163715.534591	0	11.840504	2002.0	0.0	1055.504163	12.026018	0.174528
F_778757_3113679	25009	29405	7.070	7.645	8.000	4.0	4.133	0.0	0.0	0.000	...	72132.0	102680.276817	86602.895054	0	11.389646	2002.0	0.0	0.000000	10.537294	-2.453324
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
M_74987_889210	25003	4545	89.442	9.172	41.429	7.0	0.000	0.0	381.2	148.681	...	72981.0	89094.137931	87422.722723	0	8.566898	2000.0	0.0	0.000000	8.879863	0.089586
F_182440_2984399	25003	53960	90.803	4.407	48.250	36.0	0.000	0.0	1836.0	3245.908	...	80388.0	88131.578947	94874.911482	0	NaN	NaN	0.0	0.000000	10.725508	-1.071956
M_74001_894053	25003	4545	95.959	7.229	152.000	31.0	0.000	0.0	1090.4	3458.145	...	72981.0	89094.137931	87422.722723	-1	NaN	NaN	0.0	23.843517	8.929642	-0.356848
M_67459_886161	25003	51580	109.997	4.771	71.333	34.0	0.000	0.0	1640.0	1523.297	...	69306.0	90927.272727	84795.126706	0	7.518052	1982.0	0.0	22.909716	9.359960	-0.341437
F_172467_2984485	25003	53960	115.045	2.960	46.857	40.0	0.000	0.0	604.4	390.769	...	80388.0	88131.578947	94874.911482	-3	NaN	NaN	0.0	0.000000	11.404854	-2.046122

	fips	csd_id	ha	slope	travel	p_wet	coast_2500	river_frontage	lake_frontage	lake_importance	...	bld_area_per_ha	ls_pred	fc_pred	p_f_loss_pred_100yr	ha_f_loss_pred_100yr	tco2_loss_pred_100yr	acq_cost_pred_usd	ce_cost_pred_usd	ce_cost_$/tco2	cum_tco2
pid
F_727177_2787269	25005	56375	4.738	0.629	22.0	5.0	0.0	0.0	0.0	0.000	...	0.000000	7.230578	-2.500678	25.006778	1.184821	473.928465	6.543278e+03	2.617311e+03	5.522587e+00	4.739285e+02
F_725608_2786151	25005	56375	2.373	1.725	22.0	6.0	0.0	0.0	0.0	0.000	...	0.000000	7.399235	-2.711937	27.119369	0.643543	257.417053	3.879222e+03	1.551689e+03	6.027917e+00	7.313455e+02
F_728494_2785591	25005	56375	1.495	1.300	42.0	0.0	0.0	0.0	0.0	0.000	...	0.000000	7.509981	-2.894286	28.942861	0.432696	173.078310	2.730136e+03	1.092054e+03	6.309598e+00	9.044238e+02
F_727322_2786475	25005	56375	2.578	1.041	22.0	32.0	0.0	0.0	0.0	0.000	...	0.000000	6.835246	-1.466861	14.668611	0.378157	151.262712	2.397686e+03	9.590746e+02	6.340456e+00	1.055687e+03
F_727190_2782144	25005	56375	1.657	1.097	12.0	2.0	0.0	0.0	0.0	0.000	...	0.000000	7.896622	-3.852843	38.528434	0.638416	255.366460	4.454324e+03	1.781730e+03	6.977148e+00	1.311053e+03
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
F_689217_3055697	25017	13135	1.198	0.600	5.5	13.0	0.0	0.0	0.0	0.000	...	16535.058431	14.083366	-0.001182	0.011822	0.000142	0.056649	1.565976e+06	6.263905e+05	1.105746e+07	5.778186e+07
M_120685_935014	25011	47835	1.155	0.000	42.0	48.0	0.0	155.6	0.0	0.000	...	0.000000	17.268920	-0.028356	0.283560	0.003275	1.310048	3.650722e+07	1.460289e+07	1.114683e+07	5.778186e+07
M_40312_885061	25003	21360	6.050	3.165	86.0	10.0	0.0	0.0	307.6	61.117	...	0.000000	11.641281	-0.000022	0.000215	0.000013	0.005206	6.878593e+05	2.751437e+05	5.285255e+07	5.778186e+07
M_226411_884436	25021	50250	5.010	4.274	10.0	0.0	0.0	0.0	0.0	0.000	...	26438.522954	14.671453	-0.000390	0.003901	0.000195	0.078179	1.179150e+07	4.716598e+06	6.033049e+07	5.778186e+07
M_215342_906477	25017	35425	1.550	2.032	12.0	61.0	0.0	0.0	0.0	0.000	...	3386.451613	13.715193	-0.000046	0.000457	0.000007	0.002831	1.402053e+06	5.608211e+05	1.981098e+08	5.778186e+07

	fips	csd_id	ha	slope	travel	p_wet	coast_2500	river_frontage	lake_frontage	lake_importance	...	tco2_loss_pred_100yr	acq_cost_pred_usd	ce_cost_pred_usd	ce_cost_$/tco2	offer3	offer4	tc	tc+ce_cost_pred_usd	tc+ce_cost_$/tco2	ce_cost_pred_usd_sum_with_transaction
pid
F_354523_3074640	25011	35285	26.969	11.842	125.00	0.0	0.0	0.0	162.0	11.252	...	8050.989018	2.485530e+05	99421.213240	12.348944	68663.0740	330090.549737	45000	144421.213240	17.938319	4.168910e+10
M_183328_915554	25027	34165	29.479	3.284	64.25	31.0	0.0	360.0	890.8	306.000	...	9807.050489	3.424065e+05	136962.585638	13.965727	75053.5340	402089.070054	45000	181962.585638	18.554262	5.810776e+10
F_518568_2860349	25027	63345	28.102	4.145	29.00	39.0	0.0	0.0	0.0	0.000	...	9667.088000	3.491053e+05	139642.130123	14.445108	64769.4896	396350.608000	45000	184642.130123	19.100078	5.676355e+10
F_401884_2970006	25015	52560	15.180	9.100	125.00	0.0	0.0	0.0	0.0	0.000	...	5004.677643	1.431493e+05	57259.701524	11.441237	40682.4000	205191.783353	45000	102259.701524	20.432825	4.652155e+10
F_506066_2886601	25027	68155	71.482	6.060	54.75	22.0	0.0	0.0	0.0	0.000	...	3351.532904	1.079113e+05	43164.514806	12.879037	189656.0424	137412.849083	45000	88164.514806	26.305729	5.684580e+10
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
M_139321_889963	25015	72880	2.359	6.500	27.50	148.0	0.0	426.0	0.0	0.000	...	0.000000	1.289157e+06	515662.866323	inf	2212.7420	0.000000	45000	560662.866323	inf	4.597590e+10
M_244856_915772	25009	37490	2.013	7.229	1.00	0.0	0.0	0.0	0.0	0.000	...	0.000000	1.589957e+06	635982.937821	inf	0.0000	0.000000	45000	680982.937821	inf	4.897034e+09
M_240297_934723	25009	7420	1.074	0.864	58.00	0.0	0.0	0.0	0.0	0.000	...	0.000000	3.828205e+05	153128.200192	inf	0.0000	0.000000	45000	198128.200192	inf	4.896196e+09
F_330501_2911704	25015	62745	6.743	3.781	32.50	7.0	0.0	0.0	0.0	0.000	...	0.000000	3.442470e+05	137698.808514	inf	4517.8100	0.000000	45000	182698.808514	inf	4.597749e+10
M_259920_930818	25009	21850	1.173	1.174	20.00	32.0	2.0	0.0	0.0	0.000	...	0.000000	4.870326e+05	194813.021562	inf	0.0000	0.000000	45000	239813.021562	inf	4.497034e+09

	coef	std err	t	P>\|t\|	[0.025	0.975]
const	9.7614	0.053	182.862	0.000	9.657	9.866
slope	-0.0274	0.002	-11.177	0.000	-0.032	-0.023
travel	-0.0007	0.000	-6.493	0.000	-0.001	-0.000
p_wet	-0.0051	0.000	-24.486	0.000	-0.005	-0.005
coast_2500	0.0300	0.001	30.207	0.000	0.028	0.032
river_frontage	-3.626e-05	7.71e-05	-0.471	0.638	-0.000	0.000
lake_frontage	-0.0007	5.07e-05	-12.922	0.000	-0.001	-0.001
pop_dens_tract_2012-2016	-3.376e-05	1.98e-05	-1.704	0.088	-7.26e-05	5.08e-06
hh_inc_avg_tract_2012-2016	1.316e-05	1.87e-07	70.301	0.000	1.28e-05	1.35e-05
p_f_2000	-0.0058	0.000	-16.027	0.000	-0.007	-0.005
p_prot_2000_500	8.67e-05	0.001	0.172	0.864	-0.001	0.001
p_d_2000_200	0.0126	0.001	18.356	0.000	0.011	0.014
p_d_2000_1000	-0.0004	0.001	-0.486	0.627	-0.002	0.001
p_d_2000_5000	0.0155	0.001	16.736	0.000	0.014	0.017
p_dev_open	0.0034	0.001	6.633	0.000	0.002	0.004
p_dev_low	0.0073	0.001	9.010	0.000	0.006	0.009
p_dev_medium	0.0146	0.001	19.662	0.000	0.013	0.016
p_dev_high	0.0040	0.001	5.138	0.000	0.002	0.006
p_crops	-0.0039	0.001	-4.439	0.000	-0.006	-0.002
p_pasture	-0.0001	0.001	-0.253	0.801	-0.001	0.001
p_grassland	0.0021	0.001	2.037	0.042	8.14e-05	0.004
bld_area_per_ha	1.866e-05	6.6e-07	28.272	0.000	1.74e-05	2e-05
2001	0.0837	0.039	2.160	0.031	0.008	0.160
2002	0.0817	0.038	2.152	0.031	0.007	0.156
2003	0.2283	0.038	6.033	0.000	0.154	0.302
2004	0.4424	0.037	12.035	0.000	0.370	0.514
2005	0.6535	0.037	17.885	0.000	0.582	0.725
2006	0.6194	0.037	16.565	0.000	0.546	0.693
2007	0.6743	0.038	17.898	0.000	0.600	0.748
2008	0.5678	0.039	14.575	0.000	0.491	0.644
2009	0.5697	0.040	14.301	0.000	0.492	0.648
2010	0.5044	0.039	13.089	0.000	0.429	0.580
2011	0.5055	0.040	12.682	0.000	0.427	0.584
2012	0.5223	0.039	13.351	0.000	0.446	0.599
2013	0.5588	0.039	14.379	0.000	0.483	0.635
2014	0.5453	0.039	13.931	0.000	0.469	0.622
2015	0.6266	0.041	15.312	0.000	0.546	0.707
2016	0.6110	0.043	14.099	0.000	0.526	0.696
2017	0.7002	0.109	6.409	0.000	0.486	0.914

	coef	std err	t	P>\|t\|	[0.025	0.975]
const	9.6012	0.191	50.254	0.000	9.227	9.976
slope	-0.0174	0.020	-0.877	0.380	-0.056	0.021
travel	-0.0016	0.001	-1.762	0.078	-0.003	0.000
p_wet	-0.0138	0.001	-10.100	0.000	-0.016	-0.011
coast_2500	0.0679	0.009	7.977	0.000	0.051	0.085
river_frontage	0.0053	0.001	9.006	0.000	0.004	0.006
lake_frontage	-0.0021	0.000	-4.445	0.000	-0.003	-0.001
pop_dens_tract_2012-2016	0.0009	0.000	4.706	0.000	0.001	0.001
hh_inc_avg_tract_2012-2016	2.146e-05	1.42e-06	15.101	0.000	1.87e-05	2.42e-05
p_f_2000	-0.0146	0.003	-5.670	0.000	-0.020	-0.010
p_prot_2000_500	0.0017	0.004	0.484	0.629	-0.005	0.009
p_d_2000_200	0.0138	0.004	3.101	0.002	0.005	0.022
p_d_2000_1000	0.0258	0.006	4.535	0.000	0.015	0.037
p_d_2000_5000	-0.0066	0.006	-1.053	0.292	-0.019	0.006
p_dev_open	0.0025	0.003	0.817	0.414	-0.003	0.008
p_dev_low	-0.0023	0.005	-0.430	0.667	-0.013	0.008
p_dev_medium	0.0095	0.006	1.580	0.114	-0.002	0.021
p_dev_high	-0.0162	0.008	-1.902	0.057	-0.033	0.000
p_crops	0.0026	0.007	0.348	0.728	-0.012	0.017
p_pasture	-0.0241	0.004	-6.016	0.000	-0.032	-0.016
p_grassland	-0.0202	0.007	-2.732	0.006	-0.035	-0.006
bld_area_per_ha	5.731e-05	5.32e-06	10.776	0.000	4.69e-05	6.77e-05
slope^2	0.0009	0.000	2.055	0.040	4.25e-05	0.002
travel^2	-2.908e-06	9.37e-07	-3.104	0.002	-4.74e-06	-1.07e-06
p_wet^2	1.542e-05	4.45e-06	3.466	0.001	6.7e-06	2.41e-05
coast_2500^2	-0.0002	7.04e-05	-3.222	0.001	-0.000	-8.88e-05
river_frontage^2	-2.876e-07	5.29e-08	-5.439	0.000	-3.91e-07	-1.84e-07
lake_frontage^2	7.263e-08	2.17e-08	3.354	0.001	3.02e-08	1.15e-07
pop_dens_tract_2012-2016^2	-1.153e-08	1.01e-08	-1.140	0.254	-3.14e-08	8.29e-09
hh_inc_avg_tract_2012-2016^2	-2.954e-11	2.85e-12	-10.383	0.000	-3.51e-11	-2.4e-11
p_f_2000^2	6.417e-06	1.33e-05	0.481	0.631	-1.97e-05	3.26e-05
p_prot_2000_500^2	-2.069e-05	2.16e-05	-0.956	0.339	-6.31e-05	2.17e-05
p_d_2000_200^2	-1.46e-05	4.05e-05	-0.361	0.718	-9.39e-05	6.47e-05
p_d_2000_1000^2	4.381e-05	6.27e-05	0.699	0.485	-7.9e-05	0.000
p_d_2000_5000^2	-9.542e-05	7.49e-05	-1.274	0.203	-0.000	5.14e-05
p_dev_open^2	-3.831e-05	1.97e-05	-1.949	0.051	-7.68e-05	2.22e-07
p_dev_low^2	-1.355e-05	4.15e-05	-0.326	0.744	-9.5e-05	6.79e-05
p_dev_medium^2	-0.0001	4.96e-05	-2.144	0.032	-0.000	-9.13e-06
p_dev_high^2	0.0002	7.41e-05	2.203	0.028	1.8e-05	0.000